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Abstract 


V 


We  present  a  problem  data  set  for  stochastic  programming,  and 
associated  real  world  applications.  The  problem  descriptions  were  col¬ 
lected  from  the  literature,  with  emphasis  on  variety  of  problem  struc¬ 
ture  and  application.  Each  problem  has  a  short  description,  math¬ 
ematical  problem  statement,  and  notational  reconciliation  to  a  stan¬ 
dard  problem  format.  In  addition,  most  problems  have  one  or  more 
corresponding  data  files  in  SMPS[1]  format. 
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1  Introduction 


Stochastic  programming  has  grown  in  importance  in  recent  decades,  because 
it  allows  the  modeler  to  accurately  represent  planning  under  uncertainty. 
With  strong  interest  in  solving  such  problems  and  in  finding  more  efficient 
solution  techniques,  there  has  arisen  a  need  for  a  test  set  of  stochastic  pro¬ 
gramming  problems. 

One  of  the  most  popular  forms  of  stochastic  programming  problems  is 
the  multistage  stochastic  linear  program  (MSSLP): 


minimize  Z{x{) 
subject  to 

where 


cjx  i  +  Q2(x  i) 

AiXi  =  bi 

x\  >0,  X\  £  Rni 


Qt(xt-i) E  [Qt(xt-u  ct,  bt,  At,  Tt)] , 

{ct,bt,At,Tt> 

Qt{%t— i  ?  kti  At,  Tt)  i= 

inf  {cjxt  +  Qt+i(xt)  :  Atxt  ~bt~  Ttxt-i,xt  >  0} 

xt€  Mn* 

t  =  2,3,...  , AT  —  1, 


and 


Qn{xn-\-,  cat,  bx,  An,  Tn) 

inf  {cL%  :  A^xn  =  b^  —  T^at-i, ^  0}- 

xNeM.nN 

(1) 


Here  and  throughout  this  paper,  boldface  letters  denote  random  variables, 

and  E  [•]  denotes  the  expected  value  taken  with  respect  to  x. 

{*} 

While  MSSLPs  are  growing  in  popularity,  many  of  the  applications  are 
proprietary,  and  therefore  the  models  are  not  publicly  available.  Test  set 
collections  of  MSSLPs  exist  [8,  9].  However,  they  need  to  be  enriched  with 
newer  applications.  Also,  it  would  be  helpful  if  the  original  applications 
are  described  in  the  notation  of  the  original  model,  and  related  to  a  unified 
notation  such  as  in  (1).  In  addition,  it  would  be  helpful  to  make  the  data  for 
the  test  problems  available  in  SMPS  [1]  format,  the  (emerging)  standard  for 
specifying  input  to  software  for  MSSLPs. 

To  address  the  above  needs,  we  have  collected  a  group  of  eleven  problem 
classes  from  a  variety  of  settings.  They  are  all  MSSLPs,  but  of  various  struc¬ 
tures  and  sizes.  Some  have  only  stochastic  right  hand  sides  (bt,  t  >  1),  while 
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some  have  stochastic  data  elsewhere.  In  some  cases,  problem  instances  were 
explicitly  stated  in  the  literature.  In  other  cases,  we  created  the  problems 
based  solely  on  the  problem  description  in  the  literature,  and  in  some  cases, 
there  is  not  yet  any  sample  problem. 

For  each  model  application,  we  present  a  problem  description,  a  con¬ 
cise  problem  statement,  and,  if  available,  a  numerical  example  given  by  the 
model  authors.  We  have  attempted  to  stay  as  close  to  the  authors’  notation 
as  possible  in  these  subsections.  Additionally,  where  feasible,  we  present  a 
notational  reconciliation,  which  shows  how  to  transform  the  notation  of  the 
problem  into  form  (1). 

Each  problem  class  may  be  used  to  generate  one  or  more  instances  of 
MSSLPs.  We  have  created  15  such  instances.  The  data  for  these  15  test 
problems,  as  well  as  12  other  test  problems  that  we  did  not  create,  are  avail¬ 
able  in  SMPS  format  [1]  from  the  authors. 

Each  section  covers  a  single  problem  class.  At  the  beginning  of  each 
section,  we  give  a  citation  to  the  original  application,  a  brief  description  of 
the  problem  structure,  and  if  applicable  the  names  of  the  SMPS  files  for  the 
associated  problem  instances. 

It  is  the  intention  of  the  authors  to  update  the  classes  of  applications  and 
the  test  problem  instances  as  new  application  areas  become  prominent,  and 
to  make  the  information  that  we  present  below  for  each  application  area,  as 
well  as  SMPS  inputs  for  each  test  problem  instance  freely  available  to  the 
stochastic  programming  research  community. 

In  that  spirit,  we  encourage  colleagues  to  submit  new  problem  data  with 
an  accompanying  description.  Such  submissions  should  include  the  following: 

1.  description  of  the  application  and  problem  notation, 

2.  problem  statement,  in  the  same  notation, 

3.  numerical  example,  if  practical, 

4.  reconciliation  to  the  notation  of  (1), 

5.  data  files  in  SMPS  format  for  each  instance,  and 

6.  optimal  solutions  for  each  instance  and  example. 
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2  Airlift  operations  scheduling 

Due  to  J.L.  Midler  and  R.D.  Wollmer  [12] 

(2  stage,  mixed  integer  linear  stochastic  problem) 

f/AIRL.sto. first 
/airlift/AIRL. cor,  /AIRL.tim,  < 

/AIRL.sto. second 


2.1  Description 

In  scheduling  monthly  airlift  operations,  demands  for  specific  routes  can  be 
predicted.  Actual  requirements  will  be  known  in  the  future,  and  they  may 
not  agree  with  predicted  requirements.  Recourse  actions  are  then  required 
to  meet  the  actual  requirements.  The  actual  requirements  are  expressed  in 
tons,  or  any  other  appropriate  measure,  and  they  can  be  represented  by  a 
random  variable.  Aircraft  of  several  different  types  are  available  for  service. 
Each  of  these  types  of  aircraft  has  its  own  restriction  on  number  of  flight 
hours  available  during  the  month. 

The  recourse  actions  available  include  allowing  available  flight  time  to  go 
unused,  switching  aircraft  from  one  route  to  another,  and  buying  commercial 
flights.  Each  of  these  has  its  associated  cost,  depending  on  the  type(s)  of 
aircraft  involved. 

Let  Fi  be  the  maximum  number  of  flight  hours  for  aircraft  of  type  i 
available  during  the  month,  and  let  (%•  be  the  number  of  flight  hours  required 
for  an  aircraft  of  type  i  to  complete  one  flight  of  route  j.  Then  if  Xij  is  the 
number  of  flights  originally  planned  for  route  j  using  aircraft  of  type  i,  the 
first  stage  constraint  is 

^  ^  Qijj Xjj  ^  F( ,  Mi.  (2) 

j 

When  taking  recourse  action,  we  are  under  the  constraint  that  we  cannot 
switch  away  more  flight  hours  from  aircraft  of  type  i  and  from  route  j,  than 
we  have  originally  scheduled  for  such.  This  leads  to  a  second  stage  constraint: 

^  ]  O'ijk^'ijk  —  UijXij .  Mi,Mj.  (3) 

Here,  Xijk  represents  the  increase  in  the  number  of  flights  for  route  k  flown 
by  aircraft  type  i,  because  of  being  switched  from  route  j.  Also,  is  the 
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number  of  flight  hours  required  for  aircraft  of  type  i  to  fly  route  k,  after 
having  been  switched  from  route  j.  Note  that  an  increase  of  xljk  flights  for 
route  k  results  in  the  cancellation  of 


flights  for  route  j,  since  lk  flights’  and  ‘j  flights’  are  not  necessarily  equal 
units. 

We  also  have  the  recourse  constraint  that  the  demand  for  each  route 
must  be  met.  Let  fry  be  the  carrying  capacity  (in  tons)  of  a  single  flight  of 
an  aircraft  of  type  i,  flying  route  j.  Then  the  load  originally  scheduled  to  be 
carried  on  route  j  (i.e.  the  ‘best  guess’  of  the  demand)  is 

i 

The  total  carrying  capacity  switched  away  from  route  j  in  the  recourse  action 
is 

6,  (aa)  (5) 

Conversely,  the  carrying  capacity  switched  to  route  j  is 

^  ^  bij%ikj'  (®) 

h^7 ^3 

If  we  let  yj  be  the  demand  for  route  j  which  is  contracted  commercially  in 
the  recourse,  and  yj  be  the  unused  capacity  assigned  to  route  j,  then  we 
may  combine  expressions  (4),  (5),  and  (6)  to  form  the  demand  constraint  for 
the  recourse1: 

E  -  E  %  fe)  ** + E  ^ + vt  -  y; = « 

Here,  dj  is  the  random  variable  representing  the  demand  for  route  j. 

1  We  believe  a  typographical  error  was  made  in  [12,  equation  (2.3)].  Specifically, 
‘  hjXijk  -  J//’  should  read  ‘  ^  bijxikj  +yf’. 

i,k^j  i,k^j 
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Finally,  let  Cy  be  the  cost  for  aircraft  type  i  to  be  initially  assigned  and 
fly  one  flight  of  route  j.  Let  Cijk  be  the  cost  for  aircraft  type  i  to  fly  one 
flight  of  route  k ,  after  having  been  initially  assigned  route  j.  Let  cj  be  the 
cost  per  ton  of  commercially  contracted  transport  on  route  j,  and  let  cj  be 
the  cost  per  ton  of  unused  capacity  on  route  j. 

2.2  Problem  statement 

The  problem  statement  combines  equations  (2),  (3),  and  (7): 


minimize  E  CijXij  +  Q{{%ij  } ) 

hJ 

subject  to  J2aijXij<Fu  Vi 

3 

%ij  ^  0,  Vi,  Vj, 

where 


subject  to 


^  ^  Q'ijk^'ijk  ^  Vi,  Vj 

k^3 

~Ybii  ( — )  xijk  +  Y  baxikj  +  yj  -yj  =  dj  -  Y  baxiv 

'  ttij  '  ij&j  i 

^ijk,  yj,  yj  >  o,  Vi,  Vj,  v/c. 

Note  that  the  variables  Xij  and  Xijk  represent  numbers  of  flights,  and 
therefore  should  be  integer  valued.  This  is  not  specified  in  the  problem 
statement,  however.  This  is  apparently  an  acceptable  compromise  to  Midler 
and  Wollmer  [12]  in  order  to  simplify  the  problem,  and  they  recommend  that 
the  user  adopt  his/her  own  rounding  scheme. 
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2.3  Numerical  example 

Midler  and  Wollmer  [12]  provide  a  small  numerical  example,  with  two  routes 
and  two  types  of  aircraft.  The  constants  are  given  as  follows: 


Carrying  capacity 

(tons) 

bn 

&21 

bi2 

b22 

50 

20 

75 

20 

Flying  hours 

Per  1 

round  trip 

an 

&21 

0-12 

®22 

24 

49 

14 

29 

Penalty  costs 

($/ton) 

C1 

C2 

Cl 

c2 

500 

250 

0 

0 

Cost  per  flight 

($) 

Cll 

021 

Cl2 

C22 

7200 

7200 

6000 

4000 

Costs  per  flight- 

switched  ($) 

Cl  12 

Cl21 

C212 

C221 

7000 

8200 

5500 

8700 

Flying  hours 

~ 

switched  flights 

®112 

«121 

Q212 

^221 

19 

29 

36 

56 

The  total  flying  time  available  is  Fx  -  F2  -  7200.  The  demand  for  route 
1,  di,  follows  a  lognormal  distribution  with  parameters  fj, i  =  1000,  cy  =  50, 
and  d2  independently  follows  a  lognormal  distribution  with  parameters  p2  = 
1500,  a2  =  300. 

The  optimal  solution  of  this  problem, 


'  Xu 

®12 

'  16.5 

23.2  ' 

.  ®21 

X22  . 

6.7 

0.0 

was  given  by  Midler  and  Wollmer  [12],  and  is  based  on  drawing  a  sample  of 
25  observations  from  each  distribution  of  dj  and  d2.  Midler  and  Wollmer 
[12]  did  not  specify  how  the  observations  were  drawn  from  the  distributions. 
Therefore,  we  were  not  able  to  exactly  replicate  this  problem.  We  have  cre¬ 
ated  two  versions:  airlift  .first  and  airlift .  second.  These  are  intended 
to  be  as  close  as  possible  to  the  original  problem  stated  here. 
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2.4  Notational  reconciliation 

To  make  this  problem  fit  the  notation  of  Problem  (1),  we  make  some  minimal 
changes.  Let  I  be  the  total  number  of  aircraft  types,  J  be  the  total  number 
of  routes,  and  set  n\  :=  (/)(J)  +  I ■  Set 


Xu 

Cll 

Xl2 

Cl2 

Xu 

C\J 

'  Fi  ' 

X21 

^21 

f2 

X\  ■= 

: 

,  Cl  := 

,  bi 

XIJ 

ClJ 

.  Fi  . 

Si 

0 

S2 

0 

SI 

0 

and 
A\  := 


Oil  •  •  •  alJ 

O21  • 

•  02  j 

0 

jlxl 

0 

0/1  • • •  O/ j 
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Note  that  the  number  of  stages,  N,  is  two.  The  recourse  vectors  are 


£112 

C112  —  Cnan2/uii 

Zll3 

C113  —  cii^m/ an 

Xu  J 

ciu  —  Ciiaiij/au 

X121 

C121  —  Ci2ai2i/ ai2 

X123 

Cl23  —  Ci2ai23/ai2 

X12  J 

C12J  —  C\2&\2jj  ^12 

XlJ(J-l) 

Six 

,  C2  := 

cij(j-i)  ~  cijaij(j-i)/au 

0 

S12 

0 

Su 

0 

vt 

4 

Vj 

4 

Vi 

ci 

.  vj  . 

cj 

and 


q/JxI 

di 


dj 


The  transition  matrix  is 


—an 

—a  12 

0 

qIJxI 

T2:= 

0 

—an 

By 

b2  •• 

■  Bj 

oJxl 

* 
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where 


Bi.= 


Hi 


b, 


i2 


The  matrix  A 2  is  defined  by 


0 

bij 


where 


'  A 

jIJxIJ 

qIJx2J 

i 

A 

B 

Q  JXIJ 

|  jJxJ  _jJxJ  j 

bJi 

1 

b]j 


u 

-  -  _l 

r  *-r 

?n 

bi2i  ' 

B  ■= 

bJii 

^122  ‘ 

'  '  ^J.72 

I— 

<-0 

- 1 

<~o 

Hj 


J 

^  ^  Q'ijk&kj 
k= 1 


and 


Here, 


_  J  {Qjikp/G'ik)  &pj  ^  j  —  ^ 


bijGjk 


A  j  ^  k 


{e^  G  R7  1  if  j  <  k 
ej-x  G  R7-1  if  j  >  k 
OgR7"1  ifj  =  A;. 
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3  Forest  planning 

Due  to  H.  Gassmann  [8] 

(Multistage,  linear  stochastic  problem) 

/stocforl  /stocforl.cor,  /stocforl .tim,  /stocforl . sto 
/stocfor2  /stocf or2 . cor,  /stocfor2.tim,  /stocf or2 . sto 
/stocfor3  /stocf or3 . cor,  /stocf or3. tim,  /stocf or3. sto 


3.1  Description 

The  job  of  a  long  range  forest  planner  is  to  decide  what  parts  of  the  forest 
will  be  harvested  when.  Important  criteria  for  such  a  decision  are  the  age  of 
the  trees,  and  the  likelihood  that  trees  left  standing  will  be  destroyed  by  fire. 

Gassmann  [8]  creates  a  set  of  K  age  classifications  of  equal  length  (e.g.  20 
years),  and  places  each  portion  of  the  forest  into  one  of  the  classes,  according 
to  the  age  of  the  trees  within.  He  also  divides  the  future  planning  horizon 
into  T  rounds,  each  with  a  time  length  equal  to  that  of  each  age  classification. 
That  is,  in  one  time  round,  any  trees  that  are  not  destroyed  or  harvested  will 
move  to  the  next  age  class. 

Let  the  vector  st  £  represent  the  total  amount  of  area  of  the  forest 
in  each  age  class  1  through  K  in  round  t ,  and  let  xt  €  R^  be  the  area  of  the 
forest  harvested  in  each  age  class  in  round  t.  Obviously,  we  cannot  harvest 
more  trees  of  any  age  than  currently  exist.  Therefore, 

xt<su  t  =  1,2, ...  ,T.  (8) 


Immediately  replanting  harvested  land  will  cause  an  area  increase  of  Qxt  in 
the  next  round,  where 


1  1  •••  1 

0  0  •••  0 


0  0  •••  0 


The  area  of  unharvested  trees  in  round  t  will  be  st  —  xt.  Of  this  area,  a 
random  proportion  pt  =  (pti,  Pt2>  •  •  •  ,  PtK)1"  G  R^  will  be  destroyed  by  fire 
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in  round  t.  Let 


Ptl  Pt2  •  •  •  PtK-1  PtK 

1  -  Pti  o  •  •  •  0  0 

Pt  =  o  1  -  Pt2  o  •  •  •  0 

0  •  •  •  0  1  -  PtK-1  1  -  PtK 

Then,  assuming  all  burned  areas  are  immediately  replanted,  and  therefore 
wind  up  in  age  class  1, 

St+i  =  Pt(st  Xt)  +  Qxt.  (9) 

The  material  balance  in  (9),  along  with  the  availability  limits  in  (8),  will  be 
constraints  in  the  problem. 

The  last  type  of  constraint  that  will  be  in  the  problem  is  of  the  form 

oiyrxt-i  <  yTxt  <  /V^t-i,  t  =  2, 3, . . .  ,T.  (10) 

Here  y  €  represents  the  yield,  and  a  and  (3  are  constants.  This  con¬ 
straint  might  represent  limits  on  how  fast  the  timber  industry  can  change  its 
purchasing  volume  from  one  time  period  to  the  next. 

The  objective  will  be  to  maximize  the  value  of  timber,  both  cut  and 
remaining  after  round  T,  subject  to  the  constraints  (8),  (9),  and  (10).  Since 
the  time  scale  of  the  problem  is  quite  large,  Gassmann  discounts  monetary 
values  in  future  round  t  to  current  monetary  scales  by  multiplying  by 
For  example,  if  each  round  is  20  years  long,  for  interest  (or  inflation)  rate  i, 
5  =  (1  —  i)20.  Therefore,  the  present  value  of  timber  harvested  in  round  t  is 

f'Vxt. 

If  v  E  is  the  value  of  the  trees  standing  after  round  T,  then  the  total 
value  of  trees  left  standing  after  round  T  and  cut  during  rounds  1  through 
Tis 

6t~1yTxt  +  6TvJsT+l. 

t= i 
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3.2  Problem  statement 

We  are  given  the  vector  e  RK ,  denoting  the  area  of  forest  covered  with 
timber  in  the  K  different  age  classes  at  the  beginning  of  time  period  1.  We 
are  also  given  y  E  M.K  the  vector  of  yields  (in  units  currency/hectacre  of 
forest  harvested),  v  e  RK  the  value  of  standing  timber  after  round  T,  the 
discount  rate  8,  and  constants  a,  f3. 

With  such  information,  the  problem  is  then  to 

maximize  yT  x  i  +  0,2{xi) 
subject  to  x\  <  Si  (11) 

x\  >  0, 

where 


Qt(xt-i) 

:=  E  {max  [5* 

{Pt-l.Pt . Pt> 

1yTxt  +  Qt+i(xt)]  ■ 

Xt 

<  St 

(12) 

St 

=  {Q  —  Pt-i^t-i  +  Pt-iSf- 

-1 

(13) 

otyTxt- 1 

<  yJxt  <  PyTxt- 1 } ,  t  =  2 

T 

J  *  *  *  7^3 

(14) 

Qt+i(xt ) 

:=  {5tutst+i  :  st+i  =  ( Q  ■ 

-  Pt)^t  +  Pt$t}  ■ 

(15) 

Equations  (11),  (12)  and  (15)  have  been  changed  slightly  from  the  prob¬ 
lem  statement  in  [8],  in  order  to  more  closely  match  the  content  of  the  ex¬ 
ample  problems  submitted  to  Netlib  by  Gassmann  [7]. 

3.3  Numerical  results 

Gassmann  [8]  reported  numerical  results  for  many  cases.  In  all  cases,  he 
assigned  the  values  shown  in  Table  1.  For  the  distribution  of  P,  Gassmann 
[8]  used  several  different  discretizations.  The  few  included  in  Table  2  are 
called  “upper  bound  discretizations”  by  Gassmann.  In  the  first  set  of  trials, 
Gassmann  [8]  found  the  constraints  to  be  too  severe.  Therefore,  he  changed 
the  problem  as  follows.  Violations  to  the  constraints 

ayTxt- 1  <  yTxt  <  pyrxt- 1 


12 


Table  1:  Values  of  parameters  used  in  Gassmann  [8] 


were  allowed,  but  penalized.  These  constraints  were  replaced  with 

ayJxt-i  -  yT xt  <  pti 
yTxt  -  (3yJ xt_i  <  pt2 
Ptl,Pt2  >  0, 


and  the  term 


T— 1 

+Pa) 

t= l 

was  added  to  the  objective.  In  all  the  numerical  results,  7  :=  50. 

Results  from  Gassmann  [8]  are  shown  in  Table  3.  Here,  a  discretization 
structure  of  i.jjj.kkk  means  that  an  i  point  discretization  was  used  in  the 
first  round,  a  j  point  discretization  was  used  in  rounds  two  through  four, 
and  a  k  point  discretization  was  used  in  rounds  five  through  seven.  The  only 
nonzero  component  of  the  optimal  xi  was  component  8. 

Problem  statements  in  MPS  format  may  be  found  at  Netlib,  at  http:// 
www.netlib.org/lp/data/  under  the  names  stocforl,  stocfor2,  and  stocfor3 

[7]- 
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Table  2:  Discretizations  used  in  Gassmann  [8] 


1  point  discretization 

Fire  Rate  0.06258 

Probability  1.000 

2  point  discretization 
Fire  Rate  0.08612 

0.04240 

Probability  0.4616 

0.5384 

3  point  discretization 
Fire  Rate  0.10499 

0.07354 

0.04240 

Probability  0.1847 

0.2769 

0.5384 

Table  3:  Results  from  Gassmann  [8] 


Discretization  Structure 

1.111.111 

1.222.222 

1.322.222  1.332.222 

1.333.222 

1.333.322 

Obj.  value 
Opt.  a;i(8) 

41,132.0 

20,495.8 

40,914.3 

20,047.9 

40,897.0  40,864.2 

20,076.9  19,952.8 

40,835.8 

19,947.4 

40,703.1 

19,726.6 

3.4  Notational  reconciliation 

To  express  the  problem  in  the  notation  of  Problem  (1),  we  define  the  slack 
variable  zt  :=  st  -  xt,  which  allows  us  to  eliminate  the  variable  st  for  t  >  1. 
The  vectors  ct  and  xt  are  defined  and  redefined,  respectively,  as 


'  -S*  1y  " 

'  xt 

0*xl 

Zt 

Ct  ■  = 

0 

,  Xt  ■= 

k 

0 

.  mt  . 

where  lt,  mt  €  R  are  slack  variables.  This  definition  for  ct  is  not  valid  for 
t  =  T,  when  we  have 


Ct  •*— 


-5T~1(y  +  6QTv)  ' 

0 

0 
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Let 


At  :=  [  IKxK  IKxK  0  0  ] 
and  bi  :=  si.  Then  for  t  =  2, 3, . . .  ,  T,  we  define 


■  jKxK 

jKxK 

0 

0  ' 

-  0KX 1  - 

yT 

0 

1 

0 

>  bt 

o 

1 

H 

0 

0 

1 

0 

and 

-Q  -Pt-i  0  0  ' 

Tt  :=  -PyJ  0  0  0. 

ayT  0  00 

Finally,  setting  N  :=  T,  we  have  expressed  the  problem  in  the  format  of 

(1). 
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4  Electrical  investment  planning 

Due  to  Louveaux  and  Smeers  [11] 

(Two-stage,  linear  stochastic  problem) 

/electric  /LandS. cor,  /Lands. tim,  /LandS. sto 


4.1  Description 

Louveaux  and  Smeers  [11]  consider  the  challenge  of  planning  investments  in 
the  electricity  generation  industry.  While  the  model  is,  in  general,  multistage, 
the  specific  example  given  in  [11]  is  two-stage.  The  general  iV-stage  stochastic 
model  will  be  developed  in  this  section  and  the  next,  with  the  specific  example 
in  the  following  section. 

In  each  stage  of  planning,  investments  in  n  different  technologies  may  be 
considered.  Technology  i  has  an  associated  random  investment  cost,  Ci,  a 
random  operating  cost,  qi,  and  an  availability  factor,  oz.  The  availability 
factor  is  the  portion  of  time  during  which  the  technology  may  be  operated. 

For  planned  capacity,  a  distinction  is  made  between  capacity  which  was 
planned  before  time  t  =  1,  and  that  which  was  planned  after  t  =  0.  (Here,  t 
is  an  integer.)  The  former,  pj,  includes  capacity  which  exists  on  the  ground  at 
the  start  of  the  simulation,  and  new  capacity  that  has  already  been  planned. 
The  latter  is  denoted  by  x{.  If  we  let  s*  be  the  total  capacity,  both  actual 
and  on  order,  planned  after  t  =  0,  then  we  have, 


and 


s;  =  s;  *  +  x\  -  x,~  , 


rt—Li 


i  =  1, . . .  ,n,  t  =  2, . . .  ,  N. 


Capacity  in  technology  i  also  has  a  construction  delay,  Az,  and  a  finite  life¬ 
time,  Li ,  from  planning  to  retirement.  The  total  capacity  for  technology  i  at 
time  t  is  then  (g\  +  s[-Ai). 

Demand  for  electricity  may  come  in  k  different  modes,  and  the  realization 
of  the  random  demand  in  each  mode  must  be  met  at  each  time  stage.  There¬ 
fore,  if  we  let  yij  be  the  production  of  electricity  in  mode  j  from  technology  i 
and  let  dj  be  the  random  demand  variable  for  mode  j  electricity,  we  require 


Ylyh  =  dl’  J  =  !>  ■  •  •  i  t  =  l,...,N. 

i= 1 
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A  production  balance  yields 
k 

^2 V%  -  ai(9i  +  s*_Ai)>  i  =  ,n. 

3= 1 

With  the  constraints  listed  so  far,  the  problem  does  not  have  relatively 
complete  recourse.  To  give  it  such,  Louveaux  and  Smeers  [11]  add  an  addi¬ 
tional  constraint2.  They  assume  there  is  a  technology,  which  is  always  called 
technology  n,  which  can  always  be  called  upon  to  meet  demand  in  an  imme¬ 
diate  way.  Therefore  An  is  always  zero.  Typically,  the  investment  cost  for 
technology  n  is  high.  To  simulate  purchased  electricity,  one  may  simply  let 
the  lifetime  Ln  =  1.  The  added  constraint  is 

k  n— 1 

a„(si  +  4)>EdS-E“<W  +  4"a').  t  =  1 . N.  (16) 

3= 1  1 

The  objective  is  to  minimize  the  expected  value  of  the  future  cost,  as 
represented  by  the  operating  and  investment  costs.  The  random  variable 

is  made  up  of  the  demands  (di,...  ,dn),  and  the  costs  (ci,...  ,cn)  and 

(qi,--  -  .Qn). 

4.2  Problem  statement 

We  have  the  following  definitions: 

n  =  number  of  available  technologies  (index  i ) 
k  =  number  of  modes  of  electricity  demand  (index  j) 

N  =  number  of  time  stages  (index  t) 

g\  =  capacity  of  i  to  exist  at  time  i,  decided  upon  before  t  =  1 
x\  =  new  capacity  of  i,  decided  at  time  t  >  0 

=  total  capacity,  both  actual  and  on  order,  planned  after  t  =  0 
c[  =  unit  investment  cost  of  i  at  time  t 

q*  =  unit  production  cost  of  i  at  time  t 

at  =  availability  factor  for  i 

Constraint  (16)  is  not  in  exactly  the  same  form  as  in  [11].  We  have  changed  it,  so  that 
x\  may  reflect  electricity  purchased  from  the  so-called  grid. 
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Li  =  life  of  i,  from  planning  to  retirement 
Aj  =  construction  time  for  i 
dj  =  electricity  demand  in  mode  j  at  time  t 
=  production  rate  from  i  for  mode  j  at  time  t 
T*  =  duration  of  mode  j  at  time  t 

£*  =  random  variable  whose  elements  are  { dj ,  Cj,  q[ ,  Vi,  j,  t}. 


We  are  given  all  elements  of  g,  T,  a,  L,  and  A,  with  An  —  0.  The  problem3 
is  to  choose  s,  x,  and  y  to 


N  n  j 

< 

k  \  - 

minimize 

E 

EE 

c\x\  + 

E«;':a 

€ 

,t= 1  i=l  ' 

3= 1  /  J 

subject  to  s] 

=  4 

4 

=  4 

-1+4- 

J 

i  =  1,  .  .  .  ,  77-, 

±vh 

i.  j  =  1, 

...  ,k, 

C-+- 

li 

J— 1 

i= 1 
k 


Y^ylj  ^  ai(9l +  4  Ai),  i  = 

i=i 

an(9‘n  +  4)  >  E  dj  -  E  “i  {Si  +  4_A‘)  •  t  = 

j= 1  i=l 

s\,x%y\j>  0,  i  =  ,n,  t=l,...,N,  j  =  l,...,k. 


(17) 


4.3  Numerical  results 

Louveaux  and  Smeers  [11]  present  a  two-stage  example,  with  k  =  3  operating 
modes  and  n  =  4  available  technologies.  Their  example  differs  from  their 
general  problem  development  in  several  ways.  There  is  no  immediate  source 
of  emergency  electricity,  as  A,  is  set  to  1  for  all  i.  Additionally,  there  is  a 
budget  constraint  of  120  in  stage  1.  Also,  c  and  q  are  not  stochastic.  All  of 
the  parameters  are  shown  in  Table  4.3. 

3The  term  x\  in  the  objective  function  in  (17)  is  written  as  s\  in  [11].  We  believe  this 
to  be  a  typographical  error  in  [11]. 
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Note  that  x2  =  0  and  y 1  =  0,  and  that  s*  =  xl.  Here  f  =  (3,  5,7)  with 
probabilities  (0.3, 0.4, 0.3),  respectively.  Also  note  that  with  all  technologies 
having  a  construction  delay  of  1,  x2  should  be  zero.  If  the  world  is  ending, 
there’s  no  reason  to  build  a  power  plant.  To  force  this  condition,  c2  may  be 
chosen  to  be  any  positive  vector.  Therefore,  the  problem  may  be  reduced  to 


minimize 

10xi  +  7x2  +  16x3  +  6x4  -T  E  [40yn  +  45y2i  T  32y3i  +  55y4i 

« 


+  24y!2  +  27y22  +  19.2y32  +  33y42 
+4yi3  +  4.5y23  +  3.2y33  +  5.5y43] 
subject  to  ELi  Vi  1  =  £ 

Eti  Va  =  3 
Eti  =  2 
Ei=l  ^1^12 

10xi  +  7x2  +  16x3  +  6x4  <  120 
x,y  >  0. 


ELi  Vij  <  *1 
E<=i  y*j  -  x* 

ELi  y-ij  < 
Ej= 1  v*i  —  Xa 


Louveaux  and  Smeers  [11]  report  the  optimal  solution  to  be 


x  = 


8  10  n 
3’4,  3  ,2 


T  T 


with  an  objective  value  of  381.853. 


4.4  Notational  reconciliation 

We  make  several  changes  to  problem  (17)  to  facilitate  its  transition  into  the 
format  of  problem  (1).  We  specify  c\ ,  q\ ,  and  d]  to  be  deterministic.  Further, 
we  force  Li  to  be  a  number  larger  than  N  for  %  ^  n,  and  Ln  :=  1.  Also,  let 


if  i^n 
if  i  =  n 
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4 

k 

m 

2 

m 

92 

cl 

[10,7,16,6]' 

c2 

luKSlHi 

S1 

ISlMliH 

1  jjTgg 

a 

[i,i,  i,i]' 

L 

[2, 2, 2, 2]'  | 

A 

[1,1, 1,1]'  ' 

d 1 

d2 

T 1 

SDBH^H 

T2 

BBB 

With  these  restrictions,  we  let 


where  z\  and  wl  are  slack  variables.  For  t  =  1,2,...  ,N,  let 


_ jnxn 

jnxn 

gnx(nfc) 

gnxn 

gnxl  " 

g  kxn 

g  kxn 

A 

gfcxn 

gfcxl 

gnxn 

0 

g(n-l)xn 

A 

glx(nfc) 

jnxn 

glxn 

gnxl  « 

-l 

where 


I”  jfcxfc  jkxk 


jkxk  j  g  j j^fcx(nfc) 


and 


A  := 


jlxn  Qlxn 
Qlxn  j^lxn 


glxn 

glxn 


glxn  glxn 


llxn 


Additionally,  for  t  =  2, 3, . . .  ,  N,  let 


' 

i 

£ 

X 

7 

o 

rnxn 

gnx(nfc) 

gnxn 

0 

e7 

j 

Tt:= 

cn  J 

gfcxn 

gfcxn 

gfcx(nfc) 

gfcxn 

0 

gnxn 

t 

gnx(nfc) 

gnxn 

0 

glxn 

T 

glx(nfc) 

glxn 

0  _ 

where 


Finally,  letting 


Q7i  X  1 

di 


di 

aidl 


5 


an9n 


EjU  dj  -  Er=i  ai9i  J 


we  have  put  the  problem  into  the  format  of  (1). 
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5  Selecting  currency  options 

Due  to  Klaassen,  Shapiro,  and  Spitz  [10] 

(Multistage,  non-stair  step,  linear  stochastic  problem) 

5.1  Description 

The  situation  described  by  Klaassen,  Shapiro,  and  Spitz  [10]  involves  a  U.  S. 
multi-national  corporation  (MNC),  which  has  significant  forecasted  revenues 
in  a  foreign  currency  (FC).  If  the  exchange  rate ,  S  (SUS/FC),  goes  down, 
the  MNC  would  face  declining  revenue  versus  the  forecast.  To  protect,  or 
hedge,  against  this  undesirable  possibility,  the  MNC  may  choose  to  purchase 
options  which  guarantee  a  certain  exchange  rate  at  some  point  in  the  future. 
The  guaranteed  exchange  rate  is  called  the  strike  price,  E. 

The  current  time  is  given  the  value  t  —  0,  and  the  time  at  which  the 
forecasted  revenue  will  be  realized  is  t  =  T.  The  amount  of  said  revenue 
is  assumed  to  be  known  with  certainty,  and  is  scaled  to  be  1  FC.  At  times 
t  =  0, 1, . . .  ,T  —  1,  the  MNC  may  decide  to  purchase  any  of  the  available 
options.  These  options  mature  at  time  t  =  T.  There  are  a  total  of  J  specific 
options  packages  available  for  purchase,  each  with  a  different  strike  price. 

Of  course,  the  exchange  rates  for  t  =  1, 2, . . .  ,  T  are  unknown  at  time  t  = 
0,  but  a  suitable  probability  distribution  can  be  constructed.  We  enumerate 
the  possible  exchange  rate  values  at  each  time  t  as  S},  Sf, . . .  ,  S^,  for  t  — 
0,1,...  ,  T.  Then  the  set  of  scenarios 

S  := 

{sequences  s  =  (So,  Si, . . .  ,  St)  :  1  <  S*  <  Nt,  Vt  =  1,2,...  ,  T}, 

is  the  set  of  all  possible  realizations  of  a  random  variable  s  =  (S0,  Si, . . .  ,  ST). 
Thus,  the  cardinality  of  S  is  YlJ-0Nt.  Each  realization  scenario  s  specifies 
the  exchange  rate  at  each  time  step,  and  has  an  associated  probability,  ps. 
For  scenario  s  €  <5,  we  denote  by  st  the  partial  realization  (So,  Si, ...  ,St). 

The  decision  to  purchase  options  at  any  time  t  will  depend  on  the  current 
exchange  rate,  and  on  the  type  and  quantity  of  options  previously  purchased. 
This,  in  turn,  depends  on  historical  exchange  rates.  Therefore,  the  decision 
variable  XStj  is  the  amount  of  option  j  purchased  at  time  t,  based  on  the 
partial  realization  st  =  (So, Si  , St).  This  purchase  costs  PStj  per  scaled 
unit  of  currency. 
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Because  decisions  to  purchase  options  are  only  available  through  time 
T  —  1,  there  is  no  decision  to  be  made  at  time  T.  This  allows  the  constraints 
associated  with  time  T  to  be  rolled  into  time  T  —  1. 

For  each  exchange  rate  scenario  s  G  S,  the  MNC  must  specify  an  ac¬ 
ceptable  effective  exchange  rate,  Qs,  which  would  include  effects  of  options 
purchased  as  well  as  the  actual  terminal  exchange  rate  St-  This  leads  to  the 
constraint 


T- 1  J 

st  +  E  E  Xsti 

t=o  j=l 


max{£j  —  St,  0}  —  (l  +  *ft5)T  *  P«tj  >  Qs,  Vs  G  S, 

(18) 


where  i^f  is  the  U.S.  interest  rate  at  time  t  for  partial  scenario  st.  The  left 
side  of  inequality  (18)  includes  the  payoff  from  options  which  are  active  at 
ST  as  well  as  the  discounted  cost  of  all  options  purchased.  Since  Qs  is  a 
function  of  random  variable  St,  it  is  random  as  well. 

Note  that  the  ST  in  (18)  implies  that  these  constraints  are  in  the  time 
stage  associated  with  time  T.  However,  there  are  no  time  T  decision  vari¬ 
ables,  so  we  can  write  these  constraints  in  the  time  T  —  1  stage.  At  that 
stage,  ST  is  “still”  random,  as  is  PSTj.  So,  the  only  stage  in  which  constraint 
(18)  occurs  is  the  stage  associated  with  time  T  —  1,  and  moreover,  there  are 
Nt  of  these  constraints. 

A  further  restriction  ensures  that  the  MNC  does  not  venture  into  the 
realm  of  foreign  currency  speculation.  That  is,  the  MNC  should  only  be  able 
to  purchase  options  to  cover  a  maximum  of  100%  of  the  forecasted  revenues. 
This  gives  the  constraint 


EE*^1  Vse<S.  (19) 

t= 0  j=l 

The  objective  of  this  exercise  is  to  minimize  the  expected  cost  of  all 
options  purchased: 


J 

"T— 1 

/  1  \tJ 

XsojPsoj  +  e 
j= i 

El 

_t=  1 

The  parameters  i^ts  must  be  estimated,  and  for  each  realization  St  of  St, 
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the  coefficient  PStj  is  calculated  by  the  following  formula,  given  in  [10]: 


Psa  (1  ci) 


r[<?ts(T-t)] 

E, 


-  (1  -  c2) 


,-[iftc(r-t)] 


(20) 


where 


MEj/St)  +  [%s  -  -  (T  ~  t)  (f ) 


ci  :=  N  { 


c2:=N 


and  Vt  is  the  volatility  of  the  exchange  rate,  as  measured  by  the  instantaneous 
standard  deviation  of  the  spot  rate  as  a  percentage  of  the  current  spot  rate. 
Here  N{x}  is  the  cumulative  standard  normal  distribution  function.  The 
foreign  interest  rate,  iftc  is  calculated  by 

■FC  S*  C1  +  _  1 

5<  '  E  [St+1|St]  ’ 

where  the  term  in  the  denominator  is  a  conditional  expected  value. 


5.2  Problem  statement 

Given  all  elements  of  ius,  V,  and  E,  and  given  a  discrete  probability  dis¬ 
tribution  for  s  and  corresponding  minimum  effective  exchange  rates  Qs,  we 
calculate  the  value  of  each  PStj  from  (20).  Then  the  problem  is  to 
minimize 


XI  XSojPSoj  +  E 

3= 1 


T~1  /  i  \tJ 

X^  (  1  _j_  JUS  )  X  XstjPstj 
Lt=i  +  /  j= i 


subject  to 


EE  XStj  <  1  Vs  E  S 


t= o  j= 1 
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T—l  J 

St  +  EE 

t= 0  j= 1 


AT^>0  Vj  =  1, . . .  ,  J;t  =  0, 1, . . .  ,T;  s  €  <S. 

5.3  Numerical  results 

Klaassen,  Shapiro  and  Spitz  [10]  present  a  four  stage  (T  =  4)  example,  with 
ius  =  0.10  and  V  =  0.11  for  all  time  periods  and  scenarios.  At  time  stages 
0, 1, 2,  and  3,  ten  different  options  are  available,  with  strike  prices 

(Ei,  E2, ...  ,  E10)  =  (0.44, 0.50, 0.57, 0.63, 0.70, 0.76, 0.83, 0.89, 0.96, 1.02). 

The  scenario  tree  for  the  exchange  rate  St  is  given  in  Figure  5.3.  Each 
branch  of  the  tree  has  equal  probability.  Therefore  we  have  the  following 
probabilities: 

Ps0  1 

ps\  1/3,  Vsi, 

Ps2  =  1/9, 
pS3  =  1/27,  Vs3, 
pSi  =  1/81,  Vs4. 

The  minimum  acceptable  effective  exchange  rates,  Qs  are  shown  in  Table 
4,  for  each  complete  scenario.  Results  are  given  in  Table  5.  Results  for 
different  values  of  Qs,  ius  and  V  are  also  given  in  [10]. 

5.4  Notational  reconciliation 

This  problem  does  not  fit  into  form  (1),  because  the  last  stage  (that  associ¬ 
ated  with  time  T  —  l)  contains  constraints  of  the  form  (18).  These  constraints 
contain  all  the  decision  variables  XStj,  not  just  XST_2j.  Therefore,  we  need 
“T-type”  matrices  to  connect  not  only  time  T  —  1  to  time  T  —  2,  but  also 
time  T  -  1  to  each  time  t  <  T  —  1.  We  denote  such  matrices  TTt- 


max{Ej  -  ST,  0}  -  (l  +  iuSts) 


T-t. 


stj 


>  Qs,  Vs  £  S 
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Table  4:  Minimum  acceptable  effective  exchange  rates 


scenario  (t  =  4) 

Target  Exchange  Rates  Qs 

1 

0.407 

2-5 

0.416 

6-15 

0.423 

16-31 

0.429 

32-50 

0.444 

51-66 

0.466 

67-76 

0.494 

77-80 

0.527 

81 

0.564 

In  fact,  the  speculation  constraint  (19)  also  contains  all  the  decision  vari¬ 
ables.  However,  with  a  trick,  we  can  make  these  constraints  fit  into  the  stair 
step  form  of  (1).  We  create  a  new  variable  Xzt  and  introduce  the  constraints 

j 

Xz0  —  ^  '  XSo j ) 
j= 1 


and 


Xzf  =  X 


z(t‘ 


•■>+£* 

j=i 


stji 


t  =  1,2,...  ,T-2. 


Constraint  (19)  then  may  be  written 

j 

xZ(t-2) + y>,r_u.  <  i. 

With  these  definitions,  we  can  define,  for  t  =  0, 1, . . .  ,  T  —  2, 


'  xatl ' 

PS(1 

Xt+ 1  := 

xseJ 

,  ct+1:=[l/(l  +  zf)r 

PstJ 

0 

(21) 
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Table  5:  Results  for  Klaassen,  Shapiro,  and  Spitz  [10] 

All  nonzero  option  purchases,  Xstj  are  shown.  Optimal  objective  value: 
0.1057 


Option  strike  prices 

Year 

Scenario 

0.44 

0.50 

0.57 

0.63 

0.70 

0.76 

0.83 

0.89 

0.96 

1.02 

2 

i 

0.03 

0.47 

2 

0.04 

3 

1 

3 

0.49 

4 

3 

0.0004 

0.09 

4 

0.06 

7 

0.21 

9 

10 

0.56 

0.03 

for  t  =  0, 1, . . .  ,  T  —  2, 

At+ 1  :=  [  -1  •  •  •  -1  1  ]  ,  bt+ 1  :=  0  e  R, 

and  for  t  =  2, 3  . . .  ,  T  —  1, 

Tt(t-i)  :==  [  0  •••0  -l]ellx(;+1). 

In  the  last  stage,  that  is  the  stage  associated  with  time  T  —  1,  we  will  use 
Xz(T-i)  as  a  slack  variable  in  (21). 

In  addition,  we  will  need  Nt  surplus  variables  yk,  for  k  =  1, 2, . . .  ,  Nr- 
That  is  because,  in  the  final  stage  we  also  have  the  Nt  constraints  (18).  For 
k  =  1, 2, . . .  ,  Nt,  define 

lsdk  :=  max  {Ej  -  STk,  0}  -  (l  +  i"3)^  PStj,  t  =  0, . . .  ,  T  -  1, 

j  1, . ..  , «7, 

where  Srk  is  realization  k  of  St,  given  a  partial  realization  St- i-  Note  that 
S Tk  is  random  (for  t  <  T  —  1),  because  it  is  dependent  on  sT_i. 

Let 


0 

0 

0 

r{t) 

7stn 

7sf21 

•  IfstJl 

0 

Tr(t+i) 

7stl2 

7s«22  • 

•  7stJ2 

0 

_  7stlNT 

'Yst'iNj' 

'  •  lstJNT 

0 
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for  t  =  0, 1, . . .  ,  T  —  2,  where 


r(t)  ■■= 


t  =  T-2 
otherwise. 


The  matrix  A y  is  defined  by 


1 

1 

1 

1 

0  •••  o' 

IfsT- 1 21 

Ifsx-lJl 

0 

At  *— 

7sr-i22 

J2 

0 

_jJxJ 

'Yst—i^Nt 

0 

and  the  right  hand  side  for  this  stage  is 


by  : — 


1 

Qsl  —  Syi 
Qs2  —  S  T2 


■QsNt  —  S  TNt 


The  decision  and  cost  vectors  for  the  final  stage  are 


ST-ll 


r  p 


ST- 1 


1 


X  ST  —  i  J 

xzT 

y\ 

2/2 


er:=  [l/(l +  «*,)] 


T- 1 


P 


ST- 1 


J 


0 

0 

0 


yj 


o 


With  all  of  the  definitions  above,  we  have  transformed  the  problem  into 
a  familiar  format.  It  is  a  non-stair  step  version  of  (MSSLP), 

minimize  Z{x\)  cjxi  +  Qzixi) 
subject  to  A\X\  —  b\ 

x\  >0,  Xi  €  Rm , 
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where 


Qt(.xt— x)  • —  m  [Qtfat— i)  Ctj  btj  At,  Tt)]  j 

{ct,bt,At,Tt} 


Qt{xt-\iCt,bt,At,Tt(t-\)) 

inf  {cjxt  +  Qt+iixt)  •  Atxt  =  bt-  Tt(t-i)xt-i,xt  >  0} 


t  =  2,3,...  ,N-1, 


Qn{%n- i)  • —  B  [Qn{%n-u  cn?  An, 

{cN’bN»AN'TNl>TN2*-*  ^NtN-l)} 

Tnx,  TN2, . . .  ,Tn(n-i))]  ) 


and 


Qn(xn-u  cn,  bN,  An,Tni,Tn2,  •  ■  ■  ,Tn(n- i)) 


N—l 


inf  {cL^at  :  A^x^  =  b^  -  Tm^t,  xt  >  0}. 
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6  Financial  planning  model 

Due  to  Carino  and  Ziemba  [4,  3] 

(Multistage,  linear  stochastic  problem) 

6.1  Description 

Carino  and  Ziemba  [4,  3]  describe  a  model  created  for  the  Yasuda  Fire  and 
Marine  Insurance  Co.,  Ltd.  (Yasuda  Kasai)  of  Tokyo  by  the  Frank  Russell 
Company  (Russell)  of  Takoma,  Washington.  The  model  is  a  comprehen¬ 
sive  investment,  liability,  and  risk  planning  tool.  It  is  a  multistage  linear 
stochastic  model  with  a  steady-state  condition  imposed  on  the  last  stage. 

The  complexity  of  the  model  is  such  that  it  cannot  be  completely  de¬ 
scribed  in  article  format.  The  model  presented  here  is  therefore  a  simplifica¬ 
tion  of  the  original  [4] ,  although  it  is  much  more  detailed  than  the  abbreviated 
model  presented  in  an  earlier  paper  [2,  Appendix]. 

Yasuda  Kasai  offers  many  types  of  insurance  policies,  which  differ  in  struc¬ 
ture  and  in  regulatory  treatment.  One  type  of  policy  is  a  traditional,  non¬ 
savings  insurance  policy.  Premiums  from  this  type  of  policy  go  to  the  Ya¬ 
suda  Kasai  general  account.  Other  policies  are  called  savings  policies.  These 
policies  are  really  two  policies  in  one,  with  part  of  the  premium  paying  for 
insurance  and  the  rest  constituting  a  deposit  for  savings.  The  insurance  por¬ 
tion  of  the  premium  goes  into  the  general  account,  and  the  rest  goes  into 
one  of  many  savings  accounts.  The  savings  accounts  are  separated  based 
on  regulations,  but  they  are  treated  the  same  in  this  model.  Therefore,  one 
savings  account  is  included  in  this  model. 

The  general  account  is  divided  into  a  general  allocatable  account  and 
a  non-allocatable  exogenous  account.  Funds  in  the  exogenous  account  may 
not  be  invested.  In  this  problem  description,  the  superscript  S  will  refer 
to  quantities  relating  to  the  savings  account,  while  G  and  E  will  refer  to 
those  relating  to  the  general  allocatable  and  exogenous  accounts,  respectively. 
Define  Vts,VtG,  and  VtE  as  the  market  value  of  the  savings,  general  and 
exogenous  accounts,  respectively,  at  the  end  of  period  t. 

Fund  allocations  are  not  only  classified  by  account,  they  are  also  classi¬ 
fied  by  investment  type  and  asset  class.  The  investment  type  indicates  how 
the  funds  are  invested.  Money  in  the  savings  and  general  accounts  may  be 
invested  either  directly  ( D )  or  indirectly.  There  are  three  possible  indirect 
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investment  types:  tokkin  funds  (T),  capital  to  foreign  subsidiaries  ( C ),  and 
loans  to  foreign  subsidiaries  ( L ).  So,  the  four  investment  types  are  D,  T ,  C 
and  L. 

In  contrast,  there  are  many  asset  classes,  such  as  domestic  bonds,  foreign 
equity,  and  real  estate.  In  theory,  each  combination  of  account,  investment 
type,  and  asset  class  may  have  its  own  fund  allocation.  However,  some  of  the 
combinations  are  prohibited  by  regulations.  All  of  the  permissible  allocations 
are  indexed,  and  Xnt  is  defined  as  the  allocation  of  funds  to  combination  n 
at  the  end  of  time  stage  t.  The  classifications  are  used  in  quite  a  flexible  way, 
so  that  n  €  loans  means  the  set  of  indexes  for  all  combinations  with  an  asset 
class  which  can  be  described  as  a  loan,  and  n  E  S  is  the  set  of  indexes  for 
all  combinations  involving  the  savings  account.  Therefore  the  market  value 
of  the  savings  account  can  be  expressed  by  the  constraint 

Vts~Y,Xnt  =  0.  (22) 

TltS 

The  market  value  of  the  general  account  is  written  similarly,  except  that 
it  includes  vf,  the  surplus  income  in  the  general  account.  The  constraint  is 
therefore 

V°~Y,X.  t-v?  =  0. 

n€G 

The  random  variables  in  this  model  have  dependence  on  various  rates  of 
return  and  other  company  projections.  They  are  defined  in  Table  6.  Each 
has  a  discrete  probability  distribution. 

The  savings  and  general  accounts  are  modeled  by  several  balance  and  flow 
equations.  For  example,  investment  income  Dt+ 1  is  defined,  for  the  savings 
account,  by 

Df+1  :=  ^nt+l-Xnt  +  (RInt+l  +  RP nt+l )Xnt  ~  IGf+1, 

neSD  n£SI 

and  for  the  general  account  by 

A+l  :=  ^  Rlnt+l^nt  +  (RInt+1  +  Rpnt+l)*nt  “  IG?+1- 
n^GD  n£GI 

Here,  SD  is  the  set  of  indices  corresponding  to  direct  type  allocations  from 
the  savings  account,  SI  is  the  set  corresponding  to  indirect  allocations  from 
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Table  6:  Random  variables  in  Russel-Yasuda  Kasai  model 


RInt+1 

income  return  of  allocation  n  from  the  end  of  t  to  the 
end  of  t  +  1 

RP  nt+l 

price  return  of  allocation  n  from  the  end  of  t  to  the  end 
of  t  T  1 

gt+x 

interest  rate  credited  to  policies  from  the  end  of  t  to  the 
end  of  t  +  1 

deposit  inflow  from  the  end  of  t  to  the  end  of  t  +  1 

principal  payments  from  the  end  of  t  to  the  end  of  t  +  1 

It+i 

interest  payments  from  the  end  of  t  to  the  end  of  t  +  1 

Lt 

total  reserve  liability  at  the  end  of  t 

Nt 

interest  portion  of  Lt 

IGt+i 

income  gap  resulting  from  the  difference  between  current 
market  yields  and  existing  loan  portfolio  cash  flows 

the  savings  account,  and  similarly  for  GD  and  GI  from  the  general  account. 
One  of  the  properties  of  the  indirect  investment  types  is  that  all  price  returns 
are  translated  into  income.  This  is  not  the  case  for  the  direct  investment  type. 
Therefore,  capital  gain 


X>p 


nt-f  1  Xnt) 


n£SD 


and 


Gf+1  :=  ]T  RPnt+lYn*, 

neGD 

includes  the  price  return  from  direct  investments. 

Let  Bf  be  the  income  accumulated  in  the  savings  account  through  the 
end  of  t,  and  wf  be  the  amount  transfered  to  the  savings  account  from  the 
general  account  at  the  end  of  t.  If  vf  is  the  amount  transfered  from  the 
savings  account  to  the  general  account  at  the  end  of  t,  then 

Bf+i  =  +  Df+i  ~  It+i  +  wt+ 1  —  uh-i- 

The  transfers,  wf  and  vf,  from  and  to  the  general  account,  are  established 
as  slack  variables  in  a  constraint.  This  constraint  expresses  the  desire  to  keep 
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accumulated  income  Bf+1  greater  than  accrued  interest  liability  Nf+1 .  Since 
Nt+i  =  Nt  +  St+iLt  -  !t+i>  the  constraint  is 

Bst  +  Df+ 1  +  wf+1  -  uf+1  =  Nts  +  gt+i^t  • 

The  excess  accumulated  income,  uf+1,  is  in  general  good,  because  it  con¬ 
tributes  to  income  before  taxes.  This  constraint  only  occurs  when  t  +  1  is  a 
fiscal  year-end  period. 

In  addition  to  the  income  constraint,  there  is  a  reserve  constraint,  which 
measures  the  total  reserve  shortfall  zf+1  or  surplus  qf+i-  These  are  established 
by 

V?  +  Gf+i  +  ®t+i  +  zt+ 1  ~  Qt+i  =  (1  +  gf+1)Lf . 

A  value  of  zf+1  >  0  represents  the  undesirable  situation  where  the  income 
cannot  meet  the  required  liability  reserve.  No  funds  need  be  transfered,  but 
a  penalty  is  assigned  in  the  objective. 

Another  constraint,  the  cash  flow  constraint,  addresses  the  unlikely  event 
that  net  pay  outs  from  the  savings  account,  Pf+1  +  lf+1  -  F®+1,  exceed  the 
market  value  of  the  savings  account  itself.  A  shortfall  yf+l  would  require  a 
transfer  from  the  general  account,  while  the  surplus  uf+1  is  a  slack  variable. 
The  constraint  is  expressed  as 

+  @t+ 1  +  Dt+i  +  wt+ 1  —  vt+ 1  +  Vt+ i  _  =  f*t+i  +  ^t+i  —  Ff+1.  (23) 

Any  surplus  from  constraint  (23)  is  equal  to  the  new  market  value  of  the 
savings  account: 

Vs  —  vs 

vt+ 1  ~  ut+ 1- 

Let  Bf+l  represent  the  income  accumulated  in  the  general  account  from 
the  beginning  of  the  fiscal  year  to  the  beginning  of  t  +  2.  Then 


0  if  t  +  1  is  a  year-end  stage, 

Bf  +  Df  otherwise, 


since  the  beginning  of  the  fiscal  year  is  at  the  beginning  of  t  +  2. 

Nonnegative  income  before  taxes  n*i  includes  any  income  from  the  gen¬ 
eral  account,  and  any  transfers  between  the  general  and  savings  accounts. 
The  calculation 


Yt+ 1  -  8t+ 1  =  Bf  +  Df+1  +  vf+1  -  wf+1 
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is  made  for  fiscal  year-end  stages  t  +  1  only.  Here,  sj+i  is  the  non-positive 
income,  should  such  a  dreadful  thing  occur. 

Of  course,  income  should  be  sufficient  to  pay  dividends  to  shareholders 
and  taxes.  To  encourage  such  outcomes,  let  rm  be  an  income  target,  vG+1 
be  the  income  in  excess  of  Ft+i,  and  wf+1  be  the  shortfall.  The  objective 
function  will  include  wG+1,  along  with  a  cost  penalty,  and  any  vf+1  >  0  will 
contribute  directly  to  net  worth  according  to  (22).  The  required  income 
constraint  is 


Ff+i  —  S(+i  +  Wt+i  ^t+i  —  rt+i. 

The  amount  wf+l  will  need  to  be  transferred  from  the  exogenous  account. 

The  net  worth  of  the  company  before  taxes  and  shareholder  dividends  is 
qf+l.  It  is  defined  by  the  constraint 

Vf  +  i  +  Df+\  +  Qt+i  ~  zt+  i  +  zt+ x  —  9t+ i  =  > 

where  zG+l  represents  a  negative  net  worth,  a  dire  situation. 

Cash  flow  must  also  be  balanced  in  the  general  account.  Taxes  are  as¬ 
sumed  to  be  a  constant  r  times  income  before  taxes.  Dividend  payments  to 
shareholders  are  included  in  P?+1-  Therefore,  we  have  the  constraint 

Vf  +  Gf+i  +  Df+1  —  rYt+ 1  +  vf+1  —  wf+1  —  yf+i  +  wf+1  + 

Vt+  i  —  ut+ i  ~  Pt+i  —  ^t+i' 

A  positive  value  for  yf+l  would  be  very  serious,  as  that  amount  would  have 
to  be  transferred  from  the  exogenous  account  to  pay  all  the  bills.  The  excess 
uf+1  is,  as  with  the  savings  account,  a  slack  variable  which  represents  the 
accumulated  market  value  of  the  general  account.  So, 

VG  -vG 

vt+ 1  —  ut+v 

The  accumulation  constraint  for  the  exogenous  account  includes  kf+1,  the 
projected  increase  in  the  exogenous  account: 

VtE+l  =  VtE -w?+l-y?+l+kf+l. 

In  addition  to  the  flow,  income  and  accumulation  constraints,  there  are 
many  other  constraints  in  the  model  by  Carino  and  Ziemba  [4].  These  con¬ 
straints  find  their  origins  in  external  and  internal  regulations  and  policies. 


36 


For  example,  since  loans  are  a  particularly  illiquid  asset  class,  an  internal 
policy  limits  the  change  in  allocations  to  loan  asset  investments  from  one 
period  to  the  next.  With  lb  and  ub  being  constants  defined  by  the  policy,  the 
constraint 


lb(  1  +  RInt+l  +  RPnt+l)-Xrct  5:  -Wit+1  5; 

ub{  1  +  RI'nt+1  +  RPnt+i)Ini,  n  £  loans 


is  added  to  the  model. 

The  objective  is  to  maximize  the  market  value  of  the  accounts  at  time 
T  and  minimize  the  costs  involved  with  shortfalls,  while  meeting  all  the 
constraints.  Let  c£t,  c^,  cf(,  cGt,  cft,  and  £t  be  the  cost  parameters  associated 
with  shortfalls  tof ,  yf,  zf,  wf,  zf,  and  yf,  respectively.  Let 


Then  the  objective  is  to  minimize  the  expected  value 


E 


-Vp  -  Vf  -  Vp  +  J](l  +  i)N{t'T)Ct  +  aCf  , 

t= 2 


where  7  is  the  discount  factor,  the  function  N(t,  T)  gives  the  number  of  years 
from  stage  t  to  stage  T,  a  is  the  discount  factor  for  the  end-effects  period, 
and  Cf  is  the  cost  of  shortfalls  for  the  end-effects  stage. 


6.2  Problem  statement 

Given  constants  lb,ub,T,kf,  and  Ft,  costs  and  given  discrete  distributions 
for  the  set  of  random  variables 

R  :=  {RInt,  RP„t,  IGst,  Fts,  Pts,  1^,  gts,  Lf,  IGg,  FG,  PtG  : 

t  =  1,2, . . .  ,T;for  all  n}, 


the  problem  is  to 
minimize  E 

R 

subject  to 


-V$  -  Vp  -  Vp  +  ]T(1  +  7 )N^T)Ct  +  aCf 

t= 2 

iy-y>„,  =  o 


(24) 


nes 
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«?  =  0 


y,°-Ex'’ 

n£G 

Df+1  =  RInt+lXnt  +  (RInt+l  +  RP nt+l)^nt  —  IG?+1 


tiGzSD 


n£SI 


Df+l  =  RInt+iXnt  +  y]  (RInt+l  +  RP nt+l)-^nt  IG 


G 

rt+l 


n£GD 


Gf+1  =  E  RP 


r*G 

W+i 


n€GJ 

nt+l-^nt 

nGSD 

=  ^2/  ^Pnt+l-^Tit 
neGD 


All  =  B?  +  All 


It+i  +  ^M-l  ^ll 


Bf  +  A+i  +  <1  -  til  =  N?  +  gts+1Lts 


P*5  +  Gf+1  +  Ali  + 


t+1 


A+i  - 


9<li 


(!  +  gt+i)Lt 


^  +  All  +  All  + 


t+1 


•  All  -  til  +  2/f+i 


yS  _  p5 

“X+l  _  rt+l 


+  If+i  -  F 


t+1 


Vs 

vt+ 1 


u 


't+ 1 


rg  — 

A+i  — 


0  if  t  +  1  is  a  year-end  stage, 

Bf  +  otherwise, 

At-i  —  st+i  —  +  A+i  +  ^li  —  f^li 

Ah  -  st+ 1  +  w!i  —  vt+i  =  A+i 

VtE  +  +  ^t+i  +  D?+i  +  9tli  —  zt+ 1  +  2tli  —  <7t+i  = 

P*G  +  +  Ali -  tAh  +  ^li  —  wt+i —  Vt+ 1  +  wtli  + 


Lf 


yf+ 1  wli —  pt+i 


G  _  t>G 


Ali 


yG  -?/G 

vt+ 1  —  “t+1 


K+1  =  K*  -  -  yf+1  +  fcf+1 

Z6(l  +  Rlnt+1  +  RPnt+i)Xnt  <  -At+l  5: 

1x6(1  +  RInt+i  +  RPnt+i)-Ait>  n  loans 
XnU  wf,  vfzf,  yf,  Vt,  «f ,  A  «t,  Al  v?,z?,q?,y? ,uf  >  0. 


(25) 

(26) 

(27) 


(28) 

(29) 


(30) 
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6.3  Numerical  results 


The  model  described  here  is  too  complex  for  us  to  create  empirical  data  at 
this  time.  Further,  the  original  creators  of  the  model  [4,  3]  did  not  provide 
specific  problem  data. 


6.4  Notational  reconciliation 

In  order  to  put  this  problem  in  the  notation  of  (1),  we  make  a  few  changes 
to  the  problem: 

1.  Assume  each  period  is  a  year-end  stage.  This  assumption  is  not  nec¬ 
essary,  but  we  are  required  to  state  which  stages  are  year-end,  and 
which  are  not.  The  result  of  this  supposition  is  the  elimination  from 
the  problem  of  the  variable  Bf  and  the  equation  (27). 

2.  The  end-effects  stage  is  eliminated.  This  eliminates  the  term  aCj  from 
the  objective. 

3.  Equations  (26)  and  (29)  are  eliminated  by  substituting  V  for  u  in  equa¬ 
tions  (25)  and  (28). 

4.  The  conditions  (30)  constraining  the  loans  are  eliminated. 

Order  the  number  of  accounts  in  any  way,  and  let  M  be  the  number  of 
accounts.  That  is,  the  index  n  runs  from  1  to  M.  Define  the  vectors 


and 
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Let  A  6  RMxM  be  the  diagonal  matrix  defined  for  sets 
*e{S,G,D,I}  by 


1  if  account  j  is  in  set  $ 
0  otherwise. 


Then  we  may  express  the  sums  from  the  problem  statement  in  Section 
6.2  in  matrix  notation.  For  example, 

RI„,+1X„,  =  (RIt+1)TAeA2>X,. 

ntGD 

To  begin  putting  the  problem  into  the  notation  of  (1),  set 


and 


(-11  xmaS)  i  o  ' 
(_lixMAG)  oi' 
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Then  for  t  =  2, 3, . . .  ,  T,  define 


Xt  1 

' 

qMxI 

Bst 

0 

Df 

0 

Gf 

0 

v? 

—StT 

flf 

0 

v? 

0 

wf 

(1 

-  l)T~tclt 

vf 

(1 

-  7)t-*4 

(1 

-  7)t-*4 

D? 

,  and  c4  := 

0 

G? 

0 

vtG 

—StT 

« i 3 

0 

0 

(1 

o 

1 

c^ 

1 

yf 

(1 

-  if-^yt 

2p 

Zt 

(1 

-  7  )T-*4 

Vp 

—StT 

Yt 

0 

.  s<  j 

0 

where 


if  t  =  T 
otherwise. 


The  remaining  assignments  necessary  are 


A°]. 


T,:=[T»  If], 
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and 


where  Af,  Af ,  Tf ,  and  TtG  are  defined  in  Figures  1,  2,  3  and  4,  respectively. 
Figure  1:  Array  Af  for  Russell- Yasuda  Kasai  example 
-  xM)asj  0001  0  0  0  0  O' 

(_l(i  xM)AG)  0000  0  0  0  0  0 

0  010000000 

0  000000000 

0  001000000 

0  000000000 

0  1  (-1)  0  0  0  1  (-1)  0  0 

Af  :=  0  0  1  0  0  0  0  (-1)  1  0 

0  Olio  (-1)  0  0  0  1 

0  011  (-1)  0  (-1)  1  1  o 

0  0  0  0  0  0  (-1)  1  0  0 

0  000000000 

0  00001000  (-1) 

0  0  0  0  0  0  1  (-1)  (-1)  0 

0  ooooooooo_ 


42 


Figure  2:  Array  Af  for  Russell- Yasud a  Kasai  example 

'000  0  000000  0 

0  0  1  (-1)  0  0  0  0  0  0  0 

000  0  000000  0 

100  0  000000  0 

000  0  000000  0 

010  0  000000  0 

000  0  000000  0 

Af  :=  0  0  0  0  0  0  0  0  0  0  0 

000  0  000000  0 

000  0  000000  0 

(-1)  00  0  00000  1  (-1) 

0  0  0  0  (-1)  1  0  0  0  1  (-1) 

1  1  0  (-1)  0  0  0  1  0  0  0 

1  1  (-1)  0  0  1  1  0  0  (-t)  0 

000  0  011010  0 
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Figure  3:  Array  T?  for  Russell- Yasuda  Kasai  example 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0  ' 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

[-RItTAsAD 
-(RIt  +  RPt)TA5A/] 

0 

0 

0 

0 

0 

0 

0 

0 

0 

[-RItTAGA° 

— (RIt  +  RPt)TAGA7] 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-RPtTA5AD 

0 

0 

0 

0 

0 

0 

0 

0 

0 

— RPtTAGA° 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(-1) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(-1) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  . 
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Figure  4:  Array  T[J  for  Russell- Yasuda  Kasai  example 


'  0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 


000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

000000  0  0 

1  0  0  0  0  0  1  0 

1  0  0  0  0  0  0  0 

0  0  0  0  0  0  (-1)  0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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7  Design  of  batch  chemical  plants 

Due  to  Subrahmanyam,  Pekny,  and  Reklaitis  [17] 
(Multistage,  mixed  integer  linear  stochastic  problem) 
/chem/chem .  cor,/ chem .  tim,chem .  sto 


7.1  Description 

Subrahmanyam,  Pekny,  and  Reklaitis  [17]  describe  the  design  of  a  batch  type 
chemical  plant  to  produce  products  for  which  we  do  not  know  the  future 
demand.  We  present  here  only  half  of  the  problem  given  in  [17],  the  “Design 
SuperProblem.” 

We  must  decide  how  many  plants  to  build,  of  what  type,  when  to  build 
them,  and  how  to  operate  them.  Therefore  the  problem  has  some  integer 
decision  variables.  Let  njt  be  the  number  of  new  units  of  equipment  type 
j  to  come  online  in  time  stage  t.  This  must  take  only  integer  values.  The 
cumulative  number  of  units  of  type  j  at  time  t  is  defined  as 

t 

Njt  ■—  ^  '  TijT ,  Vj,  t. 

r=  1 

The  various  plants  can  perform  different  tasks,  which  are  indexed  as  i  — 
1, . . .  ,  I.  Certain  plants  can  perform  more  than  one  task.  Let  /]asks  be  the 
set  of  tasks  from  1  to  I  which  can  be  performed  by  a  plant  of  type  j,  and  let 
jeqmp  set  of  plant  types  from  1  to  J  which  can  perform  task  i. 

We  must  decide  which  tasks  to  perform  on  which  equipment  during  each 
time  period.  Let  be  the  number  of  times  task  i  is  performed  on  equipment 
type  j  during  time  stage  t.  If  pij  is  the  processing  time  for  task  i  with 
equipment  type  j,  and  if  Ht  is  the  length  of  stage  t,  then 

Yj  PijVijt  <  HtNjU  V),  t. 

i£jt&sks 

This  constraint  enforces  the  fact  that  time  is  limited  in  each  stage.  Something 
else  which  might  limit  the  number  of  batches  is  the  much  feared  operating 
expense  budget.  Let  C. °  be  the  total  operating  budget  per  plant  for  time  t, 
and  let  c°jt  be  the  operating  expense  incurred  for  using  equipment  type  j  to 
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perform  task  i  in  stage  t.  Then  the  operating  expense  constraint  is 

E  E  <%.»««  <cr  E^-  V(- 

j  ie/)asks  i 

The  material  balance  on  the  system  includes  inventory,  production,  con¬ 
sumption,  sales,  and  purchasing  effects.  Let  Bijt  be  the  amount  of  task  i 
performed  on  plant  type  j  in  stage  t ,  measured  in  somewhat  arbitrary  reac¬ 
tion  units.  If  fSi  is  the  stoichiometric  ratio  representing  mass  of  resource  s 
produced  per  unit  of  reaction  i,  then  the  amount  of  s  produced  in  stage  t  is 

E  E  fsiBijt- 

i  ie/eqUip 

Note  that  fsi  is  negative  if  resource  s  is  consumed  in  task  i.  The  mass  of 
resource  s  in  inventory  at  the  end  of  stage  t  is  Ast,  with  maximum  limit  A™x. 
The  material  balance  constraints  are  then 

Ast  —  1)  T  ^  ^  y  ''j  fsiBijt  —  Qst  T  Qsti  t, 

i  jg/e<iui  p 

and 

Ast<A™r,  Vs,t, 

where  qsst  is  the  mass  of  resource  s  sold  in  stage  t,  and  qbst  is  that  bought  in 
stage  t. 

The  relationship  between  Bijt  and  yijt  is 

Bijt  —  mijVijt, 

Here,  is  the  capacity  of  equipment  type  j  to  perform  task  i,  measured  in 
units  of  reaction  per  batch. 

There  are  a  couple  ways  to  limit  purchases.  One  is  to  simply  impose  a 
limit,  as  in 

Qst  <  Qbst>  Vs> 

for  some  constant  Qbst.  Another  is  to  limit  capital  expenditures  to  not  exceed 
a  constant  MCt,  as  in 

ECil%  +  E<9tSAfCl,  Vi. 

j  S 
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The  symbol  Cjt  is  the  capital  investment  cost  for  a  plant  of  type  j  in  stage 
t.  The  term  vbst  is  the  value  of  purchased  resource  s  in  stage  t  per  unit  mass. 

One  of  the  random  variables  in  this  problem  is  the  demand,  Qskt,  for 
resource  s  at  time  t.  The  index  k  determines  the  scenario.  The  other  random 
variable  is  v®kt,  the  price  per  unit  mass  of  resource  s  sold  at  or  below  the 
demand  level  Q®kt  in  stage  t.  We  define  a  new  random  variable  by  rskt  := 

(Qskt>Vskt)' 

The  recourse  variables  are  q*kt,  the  amount  of  s  sold  in  stage  t  which  does 
not  exceed  demand,  and  qsskt,  the  amount  which  exceeds  demand.  Essentially, 
qss£t  is  given  away,  rather  than  sold,  as  no  credit  toward  profit  may  be  taken 
for  this  quantity.  The  recourse  variables  are  limited  by  the  equation 

Qlt  =  Qfkt  +  Qttv  Vs>M, 


and  the  inequality 

qSskt  —  Qskt)  Vs,k,t. 

In  some  industries,  it  is  important  that  the  demand  be  met  exactly.  For 
such  cases,  define  the  variable 

x  _  1 1  if  for  each  s,  qfkt  =  Q®kt, 

|0  otherwise. 

Then  we  can  set  a  guarantee  index  Gt,  to  serve  as  the  minimum  number  of 
scenarios  for  which  the  demand  may  be  met  exactly.  We  get  the  constraint 

^2xkt  >  Gt,  Vt. 

k 


The  objective  function, 


E 

t=i 


E 


s  J  /  i 

5^(vskt&1fct  -  vstQst)  -  ni(t+s)Gjt  +  '^2c°ijtyijt 

j= 1  V  i=l 


,  S  — 1 


is  the  net  present  value  of  the  facilities,  and  includes  potential  income,  capital 
expenditures  and  operating  expense.  It  contains  no  first  stage  objective 
terms,  and  should  obviously  be  maximized.  Note  that  c°Jt  =  0  if  we  cannot 
perform  task  i  with  equipment  j.  Here,  S  is  the  construction  delay  once  the 
plant  has  been  ordered. 
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7.2  Problem  statement 

We  are  given  a  discrete  probability  distribution 


{{Pkt, r.kt)  :  k  =  1,2,...  ,K}. 

In  addition,  we  require  constants  c?t,  vbst,  Cjt ,  Pij,  Ht,  C°,  fsi,  Qbst,  A™3*,  rriij , 
Gt ,  MCt,  and  5,  and  the  index  sets  /jtasks  and  7®qmp.  Then  our  goal  is  to 
maximize 


E 

t= 1 


E 


y^(vskt Qskt  ^stQst)  (  nj(M)Cjt  +  cijtVijt 


5  =  1 


3= 1 


2  =  1 


subject  to 


^  ^  PijVijt  —  HfNjt,  Vj,  t 


i(Zj  tasks 


E  E  CijtVijt  <C?J2Nit,  Vt  (31) 

j  jg/tasks  j 


t 

Njt  =  Y^njr,  Vj,t  (32) 

T  =  1 

At  =  As(t-l)  +  E  E  fsiBijt  ~  qsst  +  VM  (33) 

*  je/jequip 

<  A3“,  Va,i 

Bijt  E  miilfcj*, 

Qst  =  Qskt  +  &  Vs>  * 

Qfkt  <  Qskt>  Vs,M 
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Qst  <  Qbst>  Vs,  t 


'£C„njt  +  '£vb„qb,t<MCl,  Vt 

j  s 

(34) 

xkt=[l  if  f°r  eaCh  S’  ^  =  Q*kt’ 
lo  otherwise 

(35) 

^2xkt>Gt,  Vt 

k 

(36) 

At,  Qst,  Qsst  >  °»  Vs>  t 

qfkt,  Qsskt  >  Vs>M 

Bijt>  0,  Vi,  j,  t 

Njt)  'IT'jti  Vijt  £  ’SA . 


7.3  Numerical  results 

Subrahmanyam,  Pekny,  and  Reklaitis  [17]  present  a  problem  with  7  =  4 
tasks,  S  —  7  resources,  7  =  3  equipment  types,  T  =  2  time  stages,  and 
K  =  2  scenarios  per  stage.  Operating  costs  are  neglected,  so  =  0  for  all 
cases,  and  (31)  is  removed  from  the  problem.  Also,  the  construction  delay  6 
and  constraints  (34),  (35),  and  (36)  are  not  included  in  the  problem. 

The  parameters  for  the  problem  are  shown  in  Tables  7,  8,  9,  and  10.  Note 
that  resources  1  and  2  must  be  purchased,  resources  4  and  7  are  sold,  and 
the  remainder  are  intermediate  resources. 

Additionally, 


A™3*  =  400  Vs,*, 


Hi  =  H2  =  80  days, 
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Table  7:  Probability  distribution  for  random  variables  (demands  not  shown 
are  zero) 


kt 

Pkt 

Q’m 

Qlkt 

Kkt 

V7  kt 

11 

0.4 

0 

0 

51 

70 

21 

0.6 

150 

200 

58 

80 

12 

0.4 

0 

0 

50 

71 

22 

0.6 

150 

200 

59 

81 

Table  8:  Parameters  for  purchased  resources 


st 

< 

* 

11 

23 

200 

12 

24 

200 

21 

25 

250 

22 

26 

250 

Pll  =  P41  =  Pl2  =  P42  =  P23  =  P33  =  4, 


jtasks  =  /tasks  = 


and 


/tasks  =  {23} 

The  optimal  objective  value  stated  in  [17]  is  3300,  with  optimal  values  of 
IV21  =  1  and  jV31  =  1. 

The  problem  chem  in  our  collection  is  an  attempt  at  recreating  this  exam¬ 
ple.  We  have  not  been  able  to  verify  that  we  have  succeeded  in  this  attempt, 
as  we  have  only  run  chem  as  a  continuous  model.  The  optimal  objective  value 
for  chem  as  a  continuous  model  is  13009.16667. 
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Table  9:  Parameters  for  equipment  types 


j 

Cj  i 

Cj  2 

rriijVi 

1 

2500 

2600 

100 

2 

3000 

3100 

200 

3 

2800 

2900 

150 

Table  10:  Stoichiometric  coefficients  fSi 


i\s 

1 

2 

CO 

4 

5 

6 

7 

1 

-1 

-1 

1 

0 

0 

0 

0 

2 

0 

0 

-1 

1 

1 

0 

0 

3 

0 

-1 

0 

0 

-1 

1 

0 

4 

0 

0 

0 

0 

0 

-1 

1 

7.4  Notational  reconciliation 


For  simplicity,  equations  (35)  and  (36)  are  removed  from  the  problem,  and  6 
is  set  to  0.  The  transition  of  notation  is  then  quite  straightforward.  Define 
the  diagonal  matrix  Glw  by 

(A*-)*:-/1 

3  I  0  otherwise. 


Then,  we  make  the  following  definitions,  and  also  introduce  slack  variables 


,u 


7. 


Plj 

.-r 

i — * 

to 

0  ' 

Pj  ■■= 

;  p  •= 

0 

Pj  _ 

5? 

_ _ 1 

Vu 

Bijt 

Vjt 

; 

]l 

: 

ii 

i 

•w 

5* 

_ 1 

1 

w- 

1 _ 

Bijt  _ 
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Then,  let 


~c°t 

o-/xl 

-Ct 

0Sxl 

QlJxl 

QSxl 


vt 

05xi 

07xl 
0 

05xl 

Q/JX1 

05xl 
o5xl 
0 

We  assign  At  according  to  Figure  5.  The  vertical  lines  separate  parts  of 
the  matrix  according  to  (37).  So,  for  example,  the  second  partition  of  the 
matrix  corresponds  will  be  multiplied  by  Nt.  All  blanks  are  zero.  Note  that 
the  double  sum  in  equation  (33)  may  be  expressed  as 

E  E  /.•*«<■ 

j  jg /tasks 

The  transition  matrix  Tt  is  defined  in  Figure  6.  There  are  only  two 
nonzero  entries.  One  is  from  (33),  and  the  other  is  from  (32),  which  may  be 
rewritten 

Njt  =  Njt- 1  +  njt. 
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Our  transition  to  the  notation  of  Problem  (1)  is  complete  if  we  set 


'  Q JXJ  ‘ 

0 

Q  JXJ 
qSxS 
^max 
:=  qijxIJ 
q  SxS 

Ql 

Qbt 

MCt 


55 


igure  5:  Array  At  fi 


8  Energy  and  environmental  planning 


Due  to  Fragniere  [5] 
(Multistage,  linear  stochastic 


problem) 

/env_det . aggr 
/env.sto.irap 


/environ/env .  cor, /env .  tim,  < 


/env. sto. loose 
/env. sto. lrge 
/env.sto.xlrge 


8.1  Description 

The  model  by  Fragniere  [5]  assists  the  Canton  of  Geneva  in  planning  its  en¬ 
ergy  supply  infrastructure  and  policies.  The  model  is  based  on  the  MARKAL  , 
(market  allocation)  model.  This  is  quite  an  extensive  model,  containing  a 
great  degree  of  realism.  Included  is  the  possibility  that  emissions  of  green¬ 
house  gases  will  be  required  to  decrease.  This  possibility  is  expressed  in  a 
discrete  random  distribution. 

The  model  includes  equilibrium  constraints,  capacity  expansion  constraints, 
demand  constraints,  production  constraints,  and  environmental  constraints. 
Energy  is  supplied  by  many  different  technologies,  including  hydro  power, 
cogeneration,  fossil  fuels,  urban  waste  incineration,  and  imported  electricity. 
Demands  are  also  classified  by  technology.  Examples  are  electricity  for  in¬ 
dustrial  use,  gas  furnaces  in  existing  houses,  and  wood  stoves  in  new  houses. 
Variables  expressed  in  upper  case  letters  are  decision  variables. 

An  energy  balance  may  be  performed  on  the  supply  grid,  for  each  energy 
type.  For  types  k  which  are  neither  electricity  nor  low  temperature  heat,  the 
balance  yields 

Y,  outki(t)Pi(t)  +  y  °utki(t)cfi(t)Ci(t) +'22iMPks(t) 

iGTCH  ieDMD  s 

ifDMD 

>  Y  inpki(t)Pi(t)  +  Y  inPki(tyfi(t)Ci{t)  +  '22EXPks(t), 

ieTCH  ieDMD  s 

i£DMD 

\/k  e  ENC,Vt  eT,  (38) 
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where  the  variables  and  sets  are  defined  in  Table  11  and  Table  12,  respec¬ 
tively.  Note  that  for  i  G  DMD ,  the  term  C'j(t)  refers  to  the  installed  delivery 
capacity,  whereas  for  production  type  technologies,  it  refers  to  the  installed 
production  capacity. 

For  electricity  and  low  temperature  (district)  heat,  the  energy  balances 
are 


V 


E  *V‘) + E  IMPELCszy{t ) 

J&ELA  s 


>  ^  ^  iftPELC,i{t)QzyPi{t) 
iePRC 


+  Y  inpELcAt)cfi{t)frj{i)zyCi(t)  +  Y  EXPELCkZy(t) 
i£DMD  k 


+  rj  Y  eiPizd(t) ,  Vz  €  Z,Vy  eY,Vt  eT, 

iGSTG 

3y~n 


and 


7  E  Mt)  >  E  WPLTH,i(i)c/i(t)Ci(t)  Y  frj(i)zy, 

ieHPL  ieDMD  yeY 

VzeZ,vteT, 


respectively. 


Table  11:  Variable  definitions 


m 

the  activity,  or  utilization,  of  technology  i ,  in  period  t 

Pizy (t ) 

the  production  of  electricity  from  technology  i,  in  pe- 

riod  t ,  season  jz,  and  part  of  the  day  y 

Piz(t) 

the  production  of  low  temperature  heat  from  technol- 

ogy  i,  in  period  t ,  season  z 

Ci(t) 

the  total  installed  capacity  of  technology  i  in  period  t 

Miz(t) 

production  lost  due  to  regular  maintenance  of  technol¬ 

ogy  i  in  season  2,  period  t 

IMPks(t ) 

imported  energy  of  type  k  in  period  t ,  from  source  s 

IMPELCszy(t) 

imported  electricity,  from  source  s,  in  period  t ,  season 

and  part  of  the  day  y 

EXPks(t) 

exported  energy  of  type  k  in  period  t,  to  destination  s 

(continued  on  the  next  page) 
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Variable  definitions  (continued) 

EXPELCszy(t ) 

exported  electricity,  to  destination  s,  in  period  t,  sea¬ 
son  z,  and  part  of  the  day  y 

outki  {t^j 

output  of  energy  type  k  in  period  t,  per  unit  activity 
from  technology  i  DMD,  or  per  unit  capacity  from 
i  G  DMD 

outik(t) 

fraction  of  demand  technology  i  which  supplies  utility 
demand  k  G  DM  in  period  t 

inpki{t) 

input  of  energy  type  k  in  period  t,  per  unit  activity 
from  technology  i  £  DMD ,  or  per  unit  capacity  from 
i  6  DMD 

cfi(t) 

mean  utilization  factor  of  the  total  installed  capacity 
for  technology  i  €  DMD  in  period  t 

Qzy 

the  fraction  of  a  year  covered  by  season  z,  part  of  the 
day  y 

j{i) 

utility  demand  category  j(i)  G  DM,  for  i  G  DMD 

frj(i)zy 

fraction  of  the  utility  demand  from  category  j  ( i )  which 
comes  in  season  z,  time  of  day  y 

the  electricity  input  required  at  night  to  produce  one 
unit  of  electricity  in  the  daytime  from  technology  i  G 
STG 

n 

efficiency  coefficient  for  electrical  distribution 

p 

efficiency  coefficient  for  low  temperature  heat  distri¬ 
bution 

7 

efficiency  coefficient  for  district  (low  temperature) 
heat  distribution 

k 

duration  of  equipment  i ,  in  time  stages 

m 

new  capacity  purchased  for  technology  i,  starting  in 
period  t 

residi(t ) 

capacity  which  existed  at  the  beginning  of  the  opti¬ 
mization  problem 

demandk(t ) 

demand  for  utility  k  G  DM  in  period  t 

afi(t ) 

availability  factor  of  technology  i  in  period  t 

fOi 

the  fraction  of  a  year  that  technology  i  is  lost  for  pro¬ 
duction,  due  to  one  unit  of  unavailability 

Ui 

conversion  factor  from  units  of  capacity  to  units  of 
production 

(continued  on  the  next  page) 
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Variable  definitions  (continued) 

er 

reserve  capacity  necessary  to  cover  daily  peak  demand 
for  electricity 

hr 

reserve  capacity  necessary  to  cover  daily  peak  demand 
for  low  temperature  heat 

pki(t) 

fraction  of  installed  capacity  for  production  technol¬ 
ogy  i,  available  to  satisfy  peak  demand  in  period  t 

epki(t ) 

fraction  of  electrical  consumption  for  production  tech¬ 
nology  i,  which  corresponds  to  peak  consumption  in 
period  t 

eUj(i)  (t) 

fraction  of  capacity  for  demand  technology  i,  which 
corresponds  to  the  peak  consumption  in  period  t 

bl 

maximum  fraction  of  nighttime  electrical  production 
from  technologies  i  €  BAS 

a 

annual  discount  rate 

n 

number  of  years  per  period 

invcosti(t ) 

cost  per  unit  investment  in  technology  i ,  period  t 

fixom,i(t) 

fixed  annual  operation  and  maintenance  costs  for  tech¬ 
nology  i,  period  t,  per  unit  capacity 

varorrii(t) 

variable  annual  operation  and  maintenance  costs,  per 
unit  production,  for  non-demand  technology  i,  period 
t 

unit  cost  of  energy  type  k,  purchased  from  source  s  in 
period  t 

costks(t ) 

cost  elc, ^ ) 

unit  cost  of  electricity,  purchased  from  source  s  in  pe¬ 
riod  t 

priceks(t) 

unit  price  of  energy  type  k,  sold  to  source  s  in  period 
t 

unit  price  of  electricity,  sold  to  source  s  in  period  t 

priceks{t ) 

co2i(t) 

carbon  dioxide  emissions  per  unit  capacity,  from  tech¬ 
nology  i,  period  t 

limit  co2{t) 

limit  imposed  on  carbon  dioxide  emissions  in  period  t 

The  capacity  of  each  technology  was  either  installed  after  the  beginning 
of  the  optimization  problem,  or  it  was  there  from  the  beginning.  From  this, 
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Table  12:  Set  definitions 


ENC 

the  energy  types,  except  electricity  ( ELC )  and  low 
temperature  heat  ( LTH ) 

T 

time  periods 

TCH 

supply  and  demand  technologies 

DMD 

demand  technologies 

DMD{k ) 

demand  technologies  which  can  only  supply  utility  de¬ 
mand  k  G  DM 

DM 

utility  demands 

Y 

parts  of  the  day  (d  for  daytime,  n  for  nighttime) 

Z 

seasons  of  the  year  (w  for  winter,  s  for  summer,  i  for 
intermediate) 

ELA 

technologies  that  produce  electricity 

PRC 

energy  production  technologies 

STG 

technologies  that  effectively  allow  the  storage  of  elec¬ 
tricity 

HPL 

technologies  which  produce  low  temperature  heat 
(LTH) 

CON 

technologies  which  produce  electricity  and/or  low  tem¬ 
perature  heat 

BAS 

electrical  production  technologies  which  produce  only 
at  a  steady  rate,  day  and  night 

C02 

technologies  which  emit  carbon  dioxide 

we  get  the  constraint 

t 

m  =  £  Ii(m)  +  residiit),  it  G  T,  Vi 

m=Max{i,t-(i+i} 

We  must  meet  the  demand  for  each  utility  in  each  round.  Thus, 


£  ci(t)  + 

ieDMD(k) 


outik(t)Ci(t)  >  demandk(t), 

ieDMD 
i£DMD(k ) 

ik  G  DM,  it  G  T. 
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Of  course,  we  cannot  produce  more  than  the  capacity.  For  general  pro¬ 
duction  technologies,  this  constraint  is 

Pi(t)<afi(t)Ci(t),  Vi  G  PRC,  Vt  €  T. 

For  technologies  that  produce  electricity,  the  production  constraint  is 

Pizy(t)  +  [ -  Jz_l  -)  <  Wzy  (1  -  [1  -  afi(t)]  fa)  Ci(t), 

\  9zd  T  Qzn  / 

Vi  G  ELA,Mz  G  Z,Vy  G  F,Vi  G  T. 

The  second  term  is  the  production  lost  due  to  maintenance. 

Similarly  for  technologies  that  produce  low  temperature  heat, 

Piz{t)  +  Miz(t)  <  Ui(qzd  +  qzn)  (1  -  [1  -  afi{t)]  /o<)  Ci{t), 

Vi  G  HPL,Mz  G  Zyt  G  T. 

The  following  constraint  pertains  to  maintenance. 

Y,Mix{t)  >  [1  -  afi{t)}[  1  -  foilmCiit),  Vi  G  CON,Vt  G  T. 

On  any  given  day,  the  peak  demand  level  is,  of  course,  higher  than  the 
daily  average  demand.  The  capacity  for  production  of  electricity  must  be 
sufficient  to  cover  peak  demands,  which  occur  during  the  day  in  both  winter 
and  summer.  The  constant  er  sets  how  much  higher  than  daily  average 
demand  levels  the  peak  can  be.  The  peak  constraint  for  electricity  is 


V 


1  +  er 


T  uipki(t)Ci(t)  +  —  Y/IMPELCszd(t ) 

Qzd 


LieELA 


> 


Y.  inpELC,i(t)epki(t)Pi(t)  +  —  Y,  EXPELCszd{t ) 
iiFZc  q*d  s 

+  ^2  inpELC,i(t)elf j(i){t)cfi(t)  (~~~) 


ieDMD 


Vz  G  {w,  s},\/t  G  T. 
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The  peak  demand  constraint  for  district  heat  is 


1  +  hr 


^2  Uipki(t)Ci(t)  > 


i€HPL 


Y  inpLTHAtmt)  (frMi>wd x /r,1‘-)  Q  W.  vteT. 

„*/-  n  n  V  Hwd-rqwn  / 


ieDMD 


where  hr  is  the  analog  to  er  for  electricity. 

Some  types  of  electrical  production  technologies,  here  called  BAS,  can 
only  operate  at  a  constant  production  level,  day  and  night.  We  may  desire 
to  limit  the  percentage  of  production  from  such  technologies,  since  they  do 
not  give  hour  to  hour  operation  flexibility.  The  upper  bound,  bl  is  used  in 
the  following  constraint: 


Efl-w+E  77/ M PELCszn (i)  -  EXPELCszn(t) 

ieBAS  s 


<  bl 


J2  P^(t)  +  Y,rlIMPELCszn{t)  -  EXP ELCszn(t) 

J&ELA  s 


Vz  G  Z,Mt  G  T. 


Fragniere  [5]  states  that  the  production  of  greenhouse  gases  is  limited, 
but  we  were  unable  to  find  an  explicitly  stated  constraint.  Therefore,  we 
propose  our  own  of  the  form 

^2  co2i(t)Ci(t)  +  6  *22  IMpELCszy{t )  <  limitco2(0>  Vi  G  T. 

ieC02  S  z€Z  yeY 

(39) 

The  second  term  on  the  left  hand  side  represents  the  possibility  of  imported 
electricity  counting  toward  the  CO2  limit.  Random  S(t)  G  (0,1)  represents 
the  probability  of  such  a  rule.  Of  course,  5(1)  =  0  with  probability  one. 

The  objective  is  to  minimize  capital  and  operating  costs,  which  can  be 
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expressed  as 


E 

teT 


1 

(1  +  a)n(t_1) 


Y  invcosti(t)Ii(t )  + 

ieTCH 


X^+qO1 : 


V  771—1 


V" _ I _ 

(1  +  a)n(t  1} 


y:  fixorrii(t)Ci(t )  +  ^  waromj(l)P,(i)+ 

_i£TCH  i£PRC 


Y  Y  var°mi  w  (0  +  E  EE  varomi(i)Pizy(t )  4- 

i£HPL  z(zZ  i£ELA  z(EZ  y€zY 


Y  Y  costks (t)IMPks (t)  +  Y  Y  X  costELC,s(t)IMPELCSZy{t)  - 

keENC  s  s  zeZyeY 


Y  Ypriceks(t)Expks(t)  - 

keENC  s 

Y  YT,PriceELC’°WEXPELCs*y^ 

s  zez  yeY 


8.2  Problem  statement 

We  present  a  problem  that  is  not  as  elaborate  as  that  created  by  Fragniere  [5]. 
It  corresponds  to  the  numerical  examples  given  in  the  “Numerical  examples” 
section. 
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Minimize 


E 


1 

(1  +Q,)n(«-1) 


invcosti(t)Ii(t )  + 

ieTCH 


^m=l 


y _ I _ 

^  (1  +  Q;)n(t-1) 


X  fixorrii(t)Ci(t)  +  X  i>aromj(t)Pi(t)+ 
JeTCH  iePRC 


EE  varomi(t)Piz(t )  +  E  EE  varorrii(t)PiZy(t)  + 

ieHPL  zez  i€ELA  zeZ  yeY 


E  E  costks(t)IMPks{t )  +  EEE  costELC,s{t)IMPELCszy{t )  — 

k€ENC  s  s  z€Z  yeY 


X  y ^priceks(t)EXPks(t )  - 

keENC  s 

E  E  Er^cMEXPEW.^t) 

s  z€Z  yeY 


subject  to 


X  outki(t)Pi(t)  +  x  outki(t)cfi(t)Ci(t) +^iMPks(t) 

ieTCH  ieDMD  s 

i<£DMD 

>  inPki{t)Pi{t)  +  X  inpki(t)cfi(t)Ci(t)  +  ^r,EXPks(t), 

ieTCH  ieDMD  S 

ifDMD 

Vk  e  ENCyt  €  T, 


V 


X  pi»v(t)  +  X  IMPELCszy(t) 

jieERA  $ 


>  X  inPELcA^zyPitt) 

iePRC 


+  X  inPELC,i(t)cfi(t)frj(i)zyCi(t )  +  X  EXPELCkzy(t) 

ieDMD  k 


+  V  X  e‘P^(*),  Vz  €  Z,Vy  eY,\/teT, 

ieSTG 

92/=n 


7  E  *.(«)  *  E  inPLTHj (t)cft(t)Ci(t)  X  frj(i)zy ) 

ieFPi  ieDMD  yeY 

VzezyteT, 
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cm  =  E  Ii(m )  +  residiit ),  it  G  T,  Vi, 

m=  Max{i,(-/i+i} 

E— - 

Ci(i)  +  2^  outik(t)Ci(t)  >  demandk(t), 


ieDMD(k)  ieDMD 

'  i^DMD(k) 


ik  e  dm, 'it  €  T, 


P<(i)  <  a/i(i)a(t),  Vi  €  PPC,  Vi  G  T, 


+  g  J  <  uiQzy  (1  -  [1  -  afi{t)}  foi)  Ci(t), 


ii  G  ELA,  iz  G  Z,iy  G  y,Vi  €  T, 


EM^)  >  [1  -  a/<(i)][l  -  Vi  €  CON, it  g  T, 


V  UiPfciWCiW  +  —  Y\lMPELC8zd{t)  > 

1  +  ef  LieSiA  . 

E  inPELC,i (i) ep/jj  (i)Pj(i)  +  —  E  ^ PELCszd (i) 


i£PRC 


+  E  mPsLC,i(^)e^/j(i) (i)c/i(i)  C'j(i), 

ieDMD  '  ^2d  V 


Vz  G  {w,  s},  Vi  G  T, 


tti:  E  uiPki{t)Ci(t)  > 

+  r  ieHPL 


E  inPLTH,i{t)cfi{t) 


ieDMD 


f^j(i)wd  “1“ 
Qwd  “I" 


CiCi),  ViGT, 
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£  *»(*) + £  rjIMPELCszn{t )  -  EXPELCszn(t) 

S 

Pizn(t)  +  r)I MP ELCszn(t )  —  EX PELCszn(t) 


i£B  AS 

<  bi 


,i£E  LA 


Vz  e  Z,Vte  T, 

co2j(t)Cj(t)  +  5(i)  J2  Yj'E, IMPELCszy{t )  <  limitCoa(*),  Vi  €  T 

ieC02  S  zeZ  y£Y 


and  for  all  t  €  T, 

Pi(t)  >  o,  V*  g  PRC ,  Ii(t)  >  0 ,  Vi  e  TCtf, 

Ci(t)  >  0,  Vz  G  TCiJ,  Pi*(i)  >  0,  Vz  G  HPL,\Jz  G  Z, 

Pizy(t)  >  0,  Vz  G  ELA,\/z  G  Z,Vy  G  F,  IMPks{t )  >  0,VA;  G  ENCVs, 

IMPELCszy(t)  >  0,  Vs,  Vz  G  Z,  Vy  G  F,  EXPks(t)  >0,Vke  ENCVs, 
EXPELCszy(t )  >  0,  Vs,  Vz  G  Z,Vy  G  F,  5(1)  =  0. 

8.3  Numerical  examples 

The  problem  created  by  Fragniere  [5]  for  the  Canton  of  Geneva  is  extremely 
large  and  complex,  and  the  input  data  format  is  not  SMPS.  Therefore,  we 
have  created  our  own  sample  problems  of  this  kind.  The  numbers  in  this 
example  are  based  on  the  authors’  judgment,  not  actual  economic  data. 

The  example  creates  a  situation  similar  to  that  experienced  in  the  United 
States,  where  oil  imports  (OIL)  are  the  largest  source  of  energy.  Other 
imports  are  coal  (COL),  natural  gas  (NGS),  propane  (PRO),  nuclear  fuel 
(NUF),  and  electricity  (ELC).  There  are  no  exports  in  this  example. 

The  energy  types  allowed  are  electricity  (ELC),  gasoline  (GAS),  coal 
(COL),  heating  oil  and  diesel  (HOL),  natural  gas  (NGS),  propane/LPG 
(LPG),  jet  fuel  (JET),  and  nuclear  fuel  (NUC).  When  one  unit  of  oil  is 
imported,  the  following  portions  of  hydrocarbon  based  energy  types  are  as¬ 
sumed  to  be  gained:  0.45  gasoline,  0.25  heating  oil/ diesel,  0.10  natural  gas, 
0.10  jet  fuel,  and  0.10  propane/LPG.  The  inequalities  of  type  (38)  must  take 
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this  into  account.  For  example,  the  inequality  balancing  natural  gas  is 

0.10/MPoil(£ )  +  IMP^q  >  wpngs.hngWc/hngOOChngCO 

+  inpN  gs.nel  (t)  c/nel  (^ )  Cnel  (f)  • 

The  available  technologies  are  listed  in  Table  16,  along  with  their  asso¬ 
ciated  coefficients  for  the  example  problem.  Other  coefficients  are  listed  in 
Table  13,  Table  14  and  Table  15. 

There  are  several  two  stage  versions  of  this  problem  in  the  test  set. 
They  differ  in  how  stochasticity  is  introduced.  The  problem  env.loose ,  us¬ 
ing  the  stochastic  file  env.sto. loose,  simply  assumes  very  non-challenging 
(i.e.  loose)  CO2  limits.  The  problem  env: aggressive  (env.sto.aggr)  sets 
aggressive  C02  limits.  Each  of  these  has  five  random  realizations,  and  the 
parameter  S(t)  takes  a  value  0  with  probability  one. 

The  problem  env:import  (env.sto. imp)  uses  the  aggressive  C02  lim¬ 
its,  and,  in  addition,  considers  the  possibility  that  imported  electricity  (IM- 
PELC)  will  be  counted  toward  such  limits  in  period  two.  That  is,  5(2)  takes 
a  nonzero  value  with  nonzero  probability.  This  problem  has  fifteen  random 
realizations. 

The  problem  env.large  (env.sto.lrge)  builds  on  env.import  by  making 
random  the  costs  of  various  energy  sources.  The  number  of  realizations  is 
8,232.  The  problem  env.xlarge  (env.sto.xlrge)  is  a  larger  version  still, 
mostly  to  test  distributed  memory  capabilities  of  the  solver. 


Table  13:  Example  problem  seasonal  coefficients 


summer  winter 


day 

night 

day 

night 

Qzy 

0.60 

0.40 

0.40 

0.60 

COSt-ELC 

5.2 

5.0 

4.8 

4.6 

f  rELC, zy 

0.35 

0.25 

0.10 

0.30 
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Table  14:  Example  problem  demands 


k 

demandkiX) 

demandk(2 ) 

ELC 

170 

230 

HHO 

30 

30 

NG 

15 

25 

GAS 

60 

80 

LPG 

3 

3 

JET 

10 

20 

Table  15:  Example  problem  coefficients 


a 

0.05 

n 

5 

V 

0.80 

6HYD 

0.10 

er 

0.20 

coston 

0.8 

COStcOAL 

0.7 

COSt^ G 

0.6 

cosipRo 

0.7 

COStjvjuF 

0.9 
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Table  16:  Example  problem  technologies  and  associated  coefficients 


OCDTfOlOHlNOiiONOO 
<3  <0  t-H  H  t-H  i — I  t-H  O  t-H  t— 4  CD 


HHlOiniOOO^lClC^NH 

oooooooooooo 


qqqqqioqqqq 


OOOLOOOOiOOOOO 

COlOWIM^oOCOM^^OOOO 


✓ - s  O  O  - - v  ^ ^  r—N  O  v  ✓~'N 

O  r- liOiOO^-vOOiOOO 
'  HHH(NOOlOH^HCO(N 

1  ^3  . 

*8  ooiotnotoooiooo 


lo  t-H  T-H  CN  OO 


'  T— I  H  CO  Tf 


qoqqqivqt^TPoqqrHoo 

cod^HddddddrHod 


oq  oo  ^  ^  io  ^ 

d>  <o>  d>  d  d> 


12  12  oq  oo 

ri  n  °  ° 


HpOOBia,  J  JffjW 


>> 

-4J  rp  _ 
O  £ 

*S  *c  .s 

t3  tj  5 
^0)  g  > 


_S  -pH 

I  1 1 

^3  O  ,  r* 

G  -O 


-g  JS  o 
«•§ 

2  S  2  , 

O  o  O  c 
^  ^ 

OJ  3  <D 

§  «3  §  . 
o  o 
A  ^3 


Oj  « 

:z!^ 
S  ft  s 

O  . — r  -+J 

i  a?  ^ 
^,3  0) 
.  <+-(  m 
o  +3  a? 
0? 


O  7j 
3  ^  (D 

s  £  s 


.  ^  II 

o  o  ^  11 

<1)  (1)  H  >W 

^  ’a)  h  m 

„o  ro 
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8.4  Notational  reconciliation 

Because  of  the  size  of  the  problem,  we  reconcile  only  the  problem  from  the 
Numerical  examples  section  to  the  format  of  (1).  In  this  section,  we  will  use 
“ELI  . . .  HYD”  to  denote  the  set  of  technologies  listed  in  Table  16,  in  the 
order  presented.  We  will  also  use  “OIL  . . .  ELC”  to  denote  the  imports, 
and  “ELC  . . .  JET”  for  the  demands  in  Table  14.  Additionally,  “WD  . . . 
SN”  will  mean  the  sequence  “WD,  WN,  SD,  SN,”  and  “CEL  . . .  HYD”  will 
stand  for  the  sequence  of  electricity  producers  “CEL,  NEL,  NUL,  HYD.” 
These  abbreviations  will  make  our  arrays  smaller  to  print. 

We  will  also  use  the  notation  e»  to  mean  the  unit  vector  in  the  ith  direction 
from  the  space  K12 . 

For  t  —  1, 2,  make  the  following  definitions: 


Jeli(^) 

r(t)invcostEu(t) 

ThydW 

r  ( t )  invcostnYv  (t) 

C'eli(^) 

s(t)fixorriEu(t ) 

Chyd  (t) 

s(t)fixomKYn(t) 

Tcel.wd  (t) 

s(t)varomcEh(t) 

P CEL,SN  (t)  j 

s(t)varomCEh(t )  _ 

-PhYD.WD  (t) 

,  Ct  := 

s(t)varorriB.Yi)(t) 

Thyd,sn(^) 

s(t)varomnYT>  (t) 

IMP0lh(t) 

s(t)costo\L 

/MPnufW 

s(t)cOStNUF 

IMPELCWD[t) 

s(t)cOStELC,WD 

IMPELCSn{t ) 

_ 

s(t)costELC,SN 
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where  r(t)  :=  (1  +  a)  n(*  x)  and  s(t)  := 
Define  the  following  matrices: 


^(1  +  c)1-)  r(t).4 


BA(t)  := 


inpo  AS,  CAR  (t)  c/cAR  (0  ^5 

^PCOL.CEL  (t)cfcEh  (0e7  +  *rcpHOL,TRK(0c/TRI<(0e10 
inpnoh,Huo(t)cfimo(t)eJ  +  *«PNGS,NEL(0c/NEL(0e8 
WPNGS.HNG  (t)  c/hnG  (0e4 
*^Plpg,hlp  (^  )  c/hLP  (0  e6 
^Pjet.air  (* )  c/air  (0  ej 

WPNUC,NUL  (^  )  c/nUL  (t) &li 


) 


jBB(t)  := 


-0.45 

0 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

-0.25 

0 

0 

0 

0 

0 

-0.10 

0 

0 

0 

0 

0 

-0.10 

0 

0 

0 

0 

0 

-0.10 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

5C(i)  :=  -v  [  /4x4  /4x4  /4x4  z4x4  ]  ,  BD(t)  :=  -17  [/4x4] 


5£(t)  := 

WPELC,ELI  (<)  c/bLI  (0  /  ^ELC/WD  WPELC.ELD  (<)  c/eLD  (0  /  ^  ELC.WD 
WPELC,ELI  (0C/eLI  (<)  /  ^ELC.SN  ^PELC.ELD  (0  c/eLD  (0  /  7~ELC,SN 


q4x  10 


> 


Bf(i)  :=  Whyd  [  04x12  Z4*4],  BH(t):= 


(d  +  e2)T 

(^3  +  eio)T 

„t 


4Due  to  an  error  of  the  author’s,  r(t)  and  s(f)  were  miscalculated  in  the  production 
of  the  SMPS  files.  The  actual  values  used  in  the  SMPS  files  were  r(l)  =  s(l)  =  1.0, 
r( 2)  =  1.914,  and  s( 2)  =  4.525. 
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BN  :=  [co2eli  ■  •  •  co2hyd]> 


BK  :=  - 


WcEL^WDeJ 

WCEL^SNCy 

^HYD^WDel2 

“HYD^SNe^ 


bk  :=  - 


<?wd/  (?wd  +  <7wn) 
<?wn/ (<?wd  +  ?wn) 
9sd/ (<7sd  +  9sn) 
9sn/ (<7sd  +  Qsn) 


J 


BL(t)  := 


n 


1  +  er 


^  [—ucELPkcEh(t)e7— 


WnELP^NEL  (^ )  ®8  —  ^NULP^NUl(^)^11  —  WHYDP^HYD(^)ei2]T  + 

/.\  I  e//ELC,WD(^)/^ELC,WD/9WD  /■„•  />\.T  , 

C/EL,(()  [  ■e//ttC3DW/rBtcWflSD  J  + 

i«PELC,ELD  (OeJ  ), 


and  the  random 


BP(i)  :=  [  6{t)  S(t)  6{t)  S(t )  ] . 

Then,  finally,  we  can  assign  At,  for  t  =  1,2,  and  T2  in  blocks.  Let 


BA(t) 

BB(t) 

BE(t) 

(BC(t)  +  BF(t)) 

BD(t) 

_j  12x12 

j  12X12 

BK 

J16X16 

BL(t) 

BM 

BN 

BP 
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and 


We  define  the  random  right  hand  side  as 


07 

04 

residEu{t) 


b(  := 


resid}iYD(t) 

—demandELc{t) 

—demanduno(t) 

—demand^Gs(t) 

—demandGAs(t) 
—demand^  pg(0 
— demand  jet  (t) 
bk 

bk 

bk 

bk 

02 

limitCo2(^) 


If  the  user  then  appends  slack  variables  in  the  blocks  corresponding  to  BA(t), 
BE(t),  BH(t),  BK(t),  BL(t )  and  BN(t),  we  will  have  the  problem  in  the 
form  of  (1). 
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9  Network  model  for  asset  or  liability  man¬ 
agement 

Due  to  J.  M.  Mulvey  and  H.  Vladimirou  [15] 

See  also  Mulvey  and  Ruszczynski  [14]. 

(Two-stage,  linear  stochastic  problem) 

{  /assets. sto. small 
/assets. sto. large 


9.1  Description 

The  management  of  assets  or  liabilities  can  be  looked  at  as  a  network  prob¬ 
lem,  where  the  asset  categories  are  represented  by  nodes,  and  transactions 
are  represented  by  arcs.  The  purchase  or  sale  of  an  asset  usually  has  fixed, 
deterministic  associated  costs,  while  the  return  on  an  investment  from  one 
stage  to  the  next  is  usually  unknown. 

Let  the  set  of  nodes  be  M,  and  let  A  be  the  set  of  arcs.  There  exists  a  set 
of  terminal  arcs  T  C  A,  over  which  the  objective  value  will  be  calculated. 
Define  the  following  notation: 

Mi=the  subset  of  arcs  associated  with  deterministic  multipliers 
and  first  stage  decisions 

,42=the  subset  of  arcs  associated  with  stochastic  multipliers  and 
first  stage  decisions 

M3=the  subset  of  arcs  associated  with  second  stage  decisions 
A/i=the  subset  of  nodes  with  deterministic  balance  equations 
A/‘2=A/'\A/i 

£)+=the  set  of  outgoing  arcs  at  node  n 
D~  =the  set  of  incoming  arcs  at  node  n 
za =flow  along  arc  a  G  A\  U  Ai 
ya=flow  along  arc  a  E  Ai 
ra=multiplier  for  arc  a  E  A 
bn=supply  or  demand  at  node  n  E  M 
Za=lower  bound  for  arc  a  G  A 
ua=upper  bound  for  arc  a  G  A. 

Then  the  problem  statement  follows  simply  from  a  material  balance  at 
each  node. 
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9.2  Problem  statement 

Given  {ra  \  a  E  Ai},  {bn  :  n  e  Afx}  and  {(la,  ua)  :  a  6  *4},  the  problem  is  to 


maximize 


'y  '  'I'a^a  *b  E  ^  ^  ra^a  -t-  y  ]  raj/a 
aeAnr  {ra,bn>  LQ€^2nr  ae^3nr 


subject  to  ^2  za  -  ^2  VaZa  n  E  M, 

aGDj  a€Dn 


E  4  E  ra,2a  X/  ra2:a  + 

aeD%n{AiUA2)  aeD~C\Ai  aeD~nA2 

J2  X  ra?/a  =  bn’  n G •A/"2’ 

aG-C^nAs  aG-D^fAs 


la  —  %a  ^  ^  ^  ^  *^-25 

la<Va<  Ua,  a  E  A3. 

9.3  Numerical  examples 

Mulvey  and  Vladimirou  [15]  did  not  provide  data  for  the  numerical  examples 
that  they  discuss  [13],  so  we  have  created  two  examples,  each  with  two  stages. 
There  are  five  nodes  in  each  stage:  checking,  savings,  certificate  of  deposit 
(CD),  cash,  and  loans,  with  initial  balances  of  100,  200,  150,  80,  and  -80, 
respectively. 

Of  course,  the  yields  are  specified  as  random.  The  smaller  problem,  using 
stochastic  file  assets. sto. small,  has  100  random  realizations,  while  the 
larger  problem,  using  assets. sto. large,  has  37,500  realizations. 

9.4  Notational  reconciliation 

Suppose  the  cardinality  of  A\  U  Ai  is  n  x,  and  that  of  A3  is  ri2 ■  Enumerate 
the  arcs  so  that  arcs  1  through  nx  are  in  ,4X  U  A2-  Reorder  the  nodes  so 
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that  the  first  Ni  are  in  the  set  M,  and  the  rest  are  in  J\f2.  For  the  set 
$  E  {Ai,A2,T,V+,V-},  define  the  diagonal  matrix  Af  E  MniXni  by 


(A?) 


1  if  a  E  $ 

0  otherwise. 


Similarly,  for  the  set  <3>  E  {As,  T,T>+,Vn},  define  the  diagonal  matrix  Af  e 
r2^2  by 


1  if  (ni  +  a)  E  $ 
0  otherwise. 


The  notation  of  (1)  requires  determinism  in  all  of  the  first  stage  coeffi¬ 
cients.  However,  since  the  only  stochastic  first  stage  coefficients  are  the  costs, 
we  can  use  the  expected  value.  For  arc  a  E  Ai  U  A2  define 

fa  :=  E  [ra]. 


Note  that  for  a  E  A\,  this  is  simply  ra.  Additionally,  define  fi  E  M”1  by 


(*i) 


ra  if  a  E  A\ 
ra  if  a  E  A2- 


Similarly,  let  r 2  E  R712  be  defined  by 

(fa)a  :=  r(ni+a)>  Va  3  (m  +  a)  E  As- 
To  fit  the  notation  of  (1),  we  make  the  following  assignments: 
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10  Cargo  network  scheduling 


Due  to  Mulvey  and  Ruszczyiiski  [14j 

(Two-stage,  mixed  integer  linear  or  nonlinear  stochastic  problem) 

( /4node_det . sto 


/cargo/4node . cor, /4node.tim,  < 


/4node . sto . 16 
/4node.sto.32 
(y4node.sto.64 


10.1  Description 

Mulvey  and  Ruszczyiiski  [14]  provide  a  two  stage  network  problem  for 
scheduling  cargo  transportation.  The  flight  schedule  is  completely  deter¬ 
mined  in  stage  one,  and  the  amounts  of  cargo  to  be  shipped  are  uncertain. 
The  recourse  actions  are  to  determine  which  cargo  to  place  on  which  flights. 
Transshipment,  getting  cargo  from  node  m  to  node  n  by  more  than  one 
flight  on  more  than  one  route,  is  allowed.  When  a  transshipment  is  made, 
cargo  must  be  unloaded  at  some  intermediate  node,  so  that  it  may  be  loaded 
onto  a  different  route  going  through  the  same  node.  Such  nodes  are  called 
transshipment  nodes.  Any  undelivered  cargo  costs  a  penalty. 

The  notation  is  introduced  in  Table  17.  A  route  tt  £  V  is  a  finite  sequence 
of  nodes  (ni,  n2, . . . ,  ni)  to  be  visited  in  the  course  of  flying  the  route. 

Table  17:  Notation 


J\f= the  set  of  nodes 
'P=the  set  of  routes 
.A=the  set  of  aircraft  types 

b(m,  n)= the  amount  of  cargo  to  be  shipped  from  node  m  to 
node  n 

c(a) =the  cost  of  an  hour  of  flight  time  for  aircraft  type  a 
h( 7r,  a)=flight  hours  required  for  aircraft  type  a  to  complete 
route  7 r 

q=the  unit  cargo  cost  for  loading  and  unloading  an 
aircraft 

p=the  unit  penalty  for  undelivered  cargo _ 

(continued  on  the  next  page) 
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Notation  (continued) 

v(jc,  j)=function  which  returns  the  jth  node  in  route  n 
/(7r)=function  which  returns  the  number  of  nodes  in  route 
7 r 

cr(n) =the  maximum  number  of  landings  allowed  in  node 

n 

d(a)= the  maximum  payload  of  an  aircraft  of  type  a 
f(m,n)= the  minimum  number  of  flights  from  node  m  to 
node  n 

/imax(a)= the  maximum  flying  hours  for  aircraft  of  type  a 
hmm(a)= the  minimum  flying  hours  for  aircraft  of  type  a 
x{i r,  a)=the  number  of  aircraft  of  type  a  assigned  to  fly  route 
7 r 

d( 7r,  m,  n)=the  amount  of  cargo  delivered  directly  from  mton 
on  route  n 

t(ir,  m,  k,  n)=the  amount  of  cargo  moving  from  mton  which  is 
moved  to  transshipment  node  k  on  route  7r 
s(Tr,k,  n)=the  amount  of  transshipment  cargo  which  is  moved 
from  transshipment  node  k  to  node  n  on  route  7T 
y(m,  n)=the  amount  of  cargo  moving  from  m  to  n  which  is 
undelivered 

z(7r,  j)=the  unused  capacity  of  leg  j  on  route  7r 
U(m,  n)={,K  <E  V  :  mn  =  v(ir,ji),n  =  v(n,j2),ji  <  h) 
Vi(n)-=={n  \  n  =  v{ir,l)} 

Vi{n)—{-n  eV:n  =  v(i r,  i(7r))} 

W(n)={n  £?:n  =  v(n  ,j)  for  some  j} 

The  first  stage  constraints  include  minimum  flight  requirements 

E  E  x(ir,a)  >  f(m,n),  Vm,n£Af, 

a€A  7 r€t/(m,n) 

and  maximum  landings  limits 

E  E  x(tt,  a)  <  cr(n),  Vn  G  M. 

a€A  7 rGW'(n) 
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Assuming  the  operation  is  cyclic,  we  must  end  the  round  in  the  same  state 
as  that  in  which  we  began  the  round.  That  is, 

Y  x^a)  =  Y  x(tt,  a),  Va  £  A,  n£j\f. 

7rGVi(n)  7r  GVj(ri) 

Flying  hours  are  limited  by 

hmm(a)  <  ^a;(7r,  a)h(7T,a)  <  hmax(a),  Va  £  A- 

For  recourse  constraints,  a  cargo  material  balance  yields 


'Y  (  d(TT,  m ,  n)  +  Yj  ^(7r>  m>  k, n)  I  +  2/(m> n)  ^  b(m,  n),  Vra,  n  £  N . 
irev  V  fceW  / 


A  balance  of  all  transshipments  which  go  through  k  and  wind  up  at  n  gives 

s(7r,  fc,  n),  Vk,n£j\f. 

n€V  m€Af  ir€V 

Finally,  consider  the  loading  and  unloading  which  must  occur  throughout 
the  course  of  a  single  route.  At  the  initial  node,  we  have 


Y  (  rf(7r>  v(7r >  !)> k)  +  s(7r>  u(7r>  !)>  k)  +  Y  ^7r’  v(7r'  1)>  k> n) 
keM  V  neSS 

=  'Y,d(a)x(ir,a)  —  z(ir,l):  V7r  £  V. 

aeA 

For  the  remaining  nodes  in  the  route,  a  payload  balance  yields 


£  d(7r,v(7r,j),k)  +  s(Tr,v(-7r,j),k)  +  Y  t(7r,  v(7r> •?')>  ^ n) 


fceAT 


-  (  d(TT,k,v(TT,j))  +  s(ir,k,v(n,j))  +  ^  t(Tr,k,v(7r,j),n) 

he M  \  neJ\f  j 

=  z(n,j  -  1)  -  z(7T,j),  W  £  V,j  =  2, . . .  (/( 7r)  -  1). 
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The  objective  is  to  minimize  the  costs  and  penalties.  Mulvey  and  Ru- 
szczyhski  [14]  provide  both  a  linear  objective  function 

minimize  Z\  =  EE  c(a)h(ir,  a)x(n,a)  + 

neVaeA 


*  «E  E 


b(m,n) 


7r£V  (m,n)67T 


d( 7r,  m,  n)  +  s(tt,  m,  n )  +  ^  t( n,  m,  n,  k) 


+P  EE  y(m,  n)  >  , 

meAf  n€Af  ) 


and  a  nonlinear  objective  function 

minimize  Z2  —  ^  c(a)h( n,  a)x(iT,  a)  + 

irev  a€A 


E  < 

b(m,n) 


*I*E  E 

ttZV  (m,n)£7r 


d(7r,  m,  n )  +  s(7r,  m,  n )  +  ^  t( tv,  m,  n,  k) 


kzN 


+P  EE  y{m,n) 

m^Af  n€Af 


where 


$(:r)  =  aexp(Px). 


(40) 


10.2  Problem  statement 

Given  $  as  either  the  identity  function  or  as  in  (40),  the  problem  is  to 
minimize  Z  =  EE  c(a)/i(7T,  o)x(tt,  a)  + 

nev  aeA 


E_, .  M*£  £ 

7r€P  (m,n)G7r 


b(m,n) 


d( 7r,  m,  n)  +  s(7r,  m ,  n)  +  ^  t( n,  m,  n,  k ) 


ke  H 


+P  EE  y(m,n )  J  >  , 

meAf  neAf  )  J 
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subject  to 


E  E  x(i r,o)>/(ra,  n),  Vm,neAf, 

a€A  7 r£{/(m,n) 

E  E  x(n,  a)  <  cr(n),  Vn  G  N, 

aeA  T[^W{n) 

y  a;(7r»a)=  y  x(i r,  a),  Vo  G  w4,  n£j\f, 

7rGVi  (n)  7r6V;(n) 

hmin(a)  <  ^2;(7r,a)/i(7r,  a)  <  /jmax(a),  Vo  G  .A, 

7T€^ 

y  I  d(7r,  m,  n)  +  y  i(7r,  m,  k,  n)  )  +  y(m,  n )  >  b(m,  n),  Vm,  neJ\f, 

nev  V  fceA/"  / 


EE«*  ,  m,  A;,  n)  =  E«('  ,  fc,  n),  \/k,n€.Af, 

ir£'P  mGA/-  ttGP 


y  [  d(ir,  v(ir,  1),  k)  +  s(tt,  u(7r,  1),  A:)  +  y  t(7r,  v(?r,  1),  A:,  n) 


n€Af 


=  y  d(a)x(ir,  a)  —  z(r,  1),  Vzr  G  P, 
aG.4 


y  (  rf(?r,t>(7r,j),fc)  +  «(7rlv(7r,i),fc)  +  y  j)>  k>n) 


keAT 


n€Af 


y  (  d(n,k,  v(ir,j))  +  s(7r,k,  v(n,j))  +  y  t(ir,  k,  v(ir,  j),  n) 
keH  V  ne M  / 

=  z(n,j  -  1)  -  z{n,j),  V7T  G  V,  j  =  2, . . .  (Z(tt)  -  1), 


a;(7r,  a),  d(7r,  m,  n),t(TT, m,  k,  n ),  s(7r,  A:,  n),  y(ro,  n),  z(7r,  j)  >  0 


x(tt,  a)  G  Z. 
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10.3  Numerical  examples 

Mulvey  and  Ruszczyriski  [14]  did  not  provide  data  for  the  numerical  examples 
that  they  discuss  [13].  Therefore,  we  have  created  some  examples  from  a  four 
node  network,  with  node  airports  A,  B,  C  and  E.  All  flights  7r  £  V  have  two 
legs.  That  is,  including  the  airport  of  origin,  there  are  three  airports  in  each 
flight.  No  direct  legs  are  allowed  between  A  and  E,  but  all  other  possibilities 
are  allowed.  Flights  are  enumerated  according  to  Table  18.  The  notation 
“ABA”  means  that  the  flight  begins  at  airport  A,  flies  to  airport  B,  and 
returns  to  airport  A. 


Table  18:  Possible  flights  7r  £  V  for  the  numerical  example 


0 

ABA 

6 

BAB 

13 

ECE 

19 

CAC 

1 

ABE 

7 

BAC 

14 

ECB 

20 

CAB 

2 

ABC 

8 

BCA 

15 

ECA 

21 

CBC 

3 

ACA 

9 

BCB 

16 

EBE 

22 

CBA 

4 

ACE 

10 

BCE 

17 

EBC 

23 

CBE 

5 

ACB 

11 

BEB 

18 

EBA 

24 

CEC 

12 

BEC 

25 

CEB 

Two  types  of  airplane  are  considered.  Type  0  plane  has  a  maximum 
payload  of  8,  maximum  flight  hours  of  480,  and  costs  5  per  flight  hour.  Type 
1  plane  has  a  maximum  payload  of  6,  maximum  flight  hours  of  240,  but  only 
costs  4  per  flight  hour.  Both  types  of  airplanes  may  have  flight  hours  as  low 
as  0.  The  unit  cost,  q,  for  loading  and  unloading  is  1.0,  and  p,  the  penalty 
for  undelivered  cargo  is  1300.  There  are  no  minimum  numbers  of  flights,  and 
the  limit  on  landings  is,  for  the  base  problem,  25  for  each  airport.  Flight 
times  for  the  two  plane  types  are  listed  in  Table  19. 
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Table  19:  Flight  times  for  numerical  example 
Airplane  Type  0  Airplane  Type  1 


A 

B 

C 

E 

A 

- 

6 

8.4 

- 

B 

6 

- 

4.8 

9.6 

C 

8.4 

4.8 

- 

6 

E 

- 

9.6 

6 

- 

A 

B 

C 

E 

A 

- 

5 

7 

- 

B 

5 

- 

4 

8 

C 

7 

4 

- 

5 

E 

- 

8 

5 

- 

10.4  Notational  reconciliation 


We  reconcile  the  problem  in  the  Numerical  examples  section  to  the  form  of 
problem  (1).  Define 


£(0, 0) 
*(1,0) 

*(25,0) 
*(0,1)  ’ 
*(1.1) 


and  hi 


h(0,i) 

MM) 


h(25,i) 


L  *(25,1)  J 


for  i  =  0, 1,  and  let  h 


hO 

hi 


We  will  make  use  of  the  incidence  matrices 


W,  Vi  E  R5x26,i  =  1, 2, 3,  defined  by 


(Vi)  mn 


1  if  node  m  is  the  ith  node  in  route  n 
0  otherwise, 


and 


H? nn 


1  if  node  m  is  in  route  n 
0  otherwise. 
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We  are  now  ready  to  define  the  stage  one  problem  parameters.  Let 


and 


w 

W 

(Vi  -  V3) 

q5x26 

hj 

q1x26 

q1  x26 

hJ2 

-hi 

q1x26 

'(N 

X 

t— H 

O 

-hi 

<r(l) 

<t(5) 

QlOXl 

hmax(  0) 
hmax(  1) 
-hmin(  0) 
-hmin(l) 


c(0)h0 

c(i)h 


Stage  two  is  a  bit  more  involved.  We  order  the  stage  two  variables  into 
x2  as  shown  below.  When  trying  to  figure  out  the  ordering  rationale  for  the 
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d,  t  and  s  variables,  it  will  help  to  look  at  Table  18. 


d(0,A,B) 

. 

d(l,A,B) 

d(2,A,B) 

t (24,  C ,  E,  B) 

t(25,C,  E,  A) 

d{25,C,  E ) 

s(0,B,A) 

s(l,B,E) 

d(0,B,A) 

d(l,B,E) 

s(25,E,B) 

d( 25,  E,  B) 

y{A,  B ) 

y(A,c) 

d(l,A,B ) 

y(A,E ) 

d{2,A,E) 

y(B,A ) 

d(4,A,E) 

d(23,C,E) 

y(E,  C) 

z(  0, 1) 

d(2h,  C,  B ) 

*(1.1) 

t(0,A,B,C ) 

t(0A,B,E) 

z(  25,1) 

t(l,A,B,C) 

t(2,A,B,E) 

A  °>  2) 

f(3,  A,  C,  B) 

j 

z(25, 2) 

The  y  variables  follow  an  ordering  we  call  “ordering  J"  the  alphabetical 
ordering  on  all  combinations  of  two  nodes.  We  will  make  use  of  the  following 
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incidence  matrices,  defined  as 


ux  e  R12x68, 

(ui)«  - 

U2  €  R12x36, 

(I4)«  := 

U3  €  R12x36, 

II 

g 

UA  €  R12X26( 

(u,h  ■= 

Us  G  R26x36, 

(Us )„  ~ 

U6  G  R26x16, 

{Ue)7Tj  •= 

and 

C/7  6  R26x36, 

ii 

g 

I 


1 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

0 


0 


if  the  ith  pair  of  J  is  served  by  the  jth 

element  of  x2, 

otherwise, 

if  the  ith  pair  of  J  is  (m,  n)  in  the  jth 

t(  7r,  m,  k,  n), 

otherwise, 

if  the  ith  pair  in  J  is  ( k,n )  in  the  jth 

t(ir,  m,  k,  n), 

otherwise, 

if  the  ith  pair  in  J  is  ( k,n )  in  the  jth 

s(tt,  k,n ), 

otherwise, 

if  the  jth  t(p,  m,  k ,  n)  has  m  =  v(-k,  1)  and 

p  =  n, 

otherwise, 

if  the  (52+ j)th  element  of  x2  is  d( n,  m,  n), 
otherwise, 


if  the  jth  t(p,m,k,n)  has  k  =  v(ir,  2)  and 

p  =  7T, 

otherwise. 


We  are  finished  putting  the  problem  into  the  form  (1)  if  we  let 


u. 

u2 

q12x26 

Jl2xl2 

q12x52 

q12x68 

U3 

-u4 

gl2xl2 

q12x52 

j26x26  q26x26  jj 

U5 

q26x26 

q26x12 

j26x26  g26x26 

_ j26x26  j26x26  q26x26 

~U7 

j26x26 

q26x12 

_j26x26  j26x26 

c2  :=  [  9(llx68)  |  ?(llx36)  |  q( llx26)  |  p(alx12)  |  0lx52  ]  ,  b2 


b 

q12x1 
q26x1  > 

q26x1 
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and 


q12x52 

q12x52 

-d(0)/26x26  -d(l)/26x26 

q26x52 
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11  Telecommunication  network  planning 

Due  to  Sen,  D overspike  and  Cosares  [16] 

(Two  stage,  mixed  integer  linear  stochastic  problem) 

/phone/phone .  cor, /phone .  tim, phone .  sto 


11.1  Description 

The  service  of  providing  private  lines  to  telecommunication  customers  is  one 
with  which  most  people  are  not  familiar.  Such  service  is  used  by  large  cor¬ 
porations  between  business  locations  for  high  speed,  private  data  transmis¬ 
sion.  Private  lines  are  generally  used  for  a  much  longer  duration  than  public 
switched  service,  and  they  generally  carry  more  capacity  per  connection. 

A  manager  of  such  a  network  must  be  constantly  looking  to  the  future,  de¬ 
ciding  where  and  how  much  to  expand  capacity.  In  this  problem  formulation, 
the  “how  much”  is  decided  beforehand,  to  some  extent,  by  the  imposition  of 
an  overall  budget.  Within  the  constraints  of  the  budget,  expansion  is  not  pe¬ 
nalized.  The  goal  is  to  minimize  the  unserved  requests,  while  staying  within 
budget. 

Such  networks  are  usually  very  interconnected,  so  that  for  any  point-to- 
point  demand  pair,  there  is  usually  more  than  one  route  which  may  service 
the  demand.  Each  route  is  made  of  one  or  more  direct  links. 

Let  n  be  the  number  of  direct  links  in  the  network  which  might  be  ex¬ 
panded,  and  let  x  E  Z"  be  the  vector  of  expanded  capacities  in  the  links, 
where  Z  is  the  set  of  integers.  Let  m  be  the  number  of  point-to-point  pairs 
to  be  served  by  the  network,  and  d  E  Zm  be  the  random  variable  of  demands 
for  service  between  the  pairs. 

The  total  budget  constraint  will  be  denoted  by  b.  Then,  the  problem  is 
to 

minimize  E  [< Q(x ,  d)] 

d 

n 

subject  to  E  xj  <  b, 

3= 1 

X  >  0, 

where  Q(x,  d)  represents  the  number  of  unserved  requests,  subject  to  network 
balance  constraints. 
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For  point-to-point  pair  i  —  1, . . .  ,m,  let  R(i)  be  the  set  of  routes  which 
may  be  used  to  satisfy  a  request  for  service  between  the  two  locations.  Ad¬ 
ditionally,  for  route  r  £  R(i),  let  air  €  Zn  be  the  incidence  vector  defined 
by 


iair)j 


1  if  link  j  £  r, 
0  otherwise. 


Let  e  £  Zn  be  the  existing  capacity  in  the  network. 

The  recourse  variables  are  s*,  the  number  of  unserved  requests,  and  fir, 
the  number  of  connections  serving  point-to-point  pair  i  over  route  r.  Then, 
the  recourse  problem  is 

m 

Q(x,d )  :=minimize  ^  Sj 

2—1 

m 

subject  to  ^2  ^22  airfir  —  x  ■**  e> 
i- 1  reR(i ) 

^  ^  fir  T  Sj  (d)j,  Vi  1,  .  .  .  ,771 
reR(i) 

fir,  Si  >  0,  Vi,reR(i) 
fir  £  Vi,  r  £  i?(i). 
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11.2  Problem  statement 


Given  the  budget  constraint  b,  and  the  current  condition  of  the  network 
{air,  e},  the  problem  is  to 


minimize 


subject  to 


E 

d 

n 

^2(x)i  < 

j=i 

m 

EE  O’irfir  E  T  C, 

i=l  r£it(i) 

^  ^  fir  T  S*  —  (^Oi)  ^  1)  •  ■  •  > 

rER(i) 

x,fir,Si>  0,  Vi,r  E  R(i) 
x,  fir  G  Z,  Vi,  r  €  .R(i). 


11.3  Numerical  example 

We  have  created  an  example  with  215  =  32, 768  random  realizations  and  six 
nodes.  The  possible  routing  is  illustrated  in  Figure  7,  and  the  possible  routes 
connecting  each  two-node  combination  are  enumerated  in  Table  20. 

The  initial  capacity  of  the  network,  e,  is  as  follows: 


route 

1 

2 

3 

4 

5 

6 

7 

8 

capacity 

2 

2 

4 

4 

2 

4 

3 

1 

11.4  Notational  reconciliation 

To  put  the  problem  into  the  notation  of  (1),  let  z\  G  E  and  z2  G  En  be  slack 
variables.  Then  set 


xi  := 


x 

Zi 


• —  n"+i 


Ci  :=  0 


Al  :=  llx(n+1),  :=  b  G  E, 
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Table  20:  Enumeration  of  all  possible  routes  for  telephone  network  example 
node  1  — >  2  node  1  — >  3  node  1  — >  4  node  1  — >  5  node  1  6 


0  12 

1  142 

2  1452 

3  1342 

4  13452 

node  2  — j  3 
0  213 

1  2143 

2  243 

3  2543 

4  25413 

5  2413 


0  13 

1  143 

2  1243 

3  12543 


node  2  — ¥  4 
0  24 

1  214 

2  2134 

3  254 


0  14 

1  124 

2  1254 

3  134 


node  2  — >  5 
0  25 

1  245 

2  2145 

3  21345 


0  125 

1  1245 

2  145 

3  1425 

4  1345 

5  13425 
node  2  — >■  6 


0  136 

1  1436 

2  12436 

3  125436 


node  3  — j  4 


0  2136 

1  21436 

2  2436 

3  24136 

4  25436 

5  254136 


0  34 

1  314 

2  3124 

3  31254 


node  3  — >  5 
0  345 

1  3425 

2  3145 

3  31425 

4  34125 

5  31245 

6  3125 


node  3  — >  6 
0  36 


node  4  — >•  5 
0  45 

1  425 

2  4125 

3  43125 


node  4  — »  6 
0  436 

1  4136 

2  42136 

3  452136 


node  5-4  6 
0  5436 

1  54136 

2  52436 

3  524136 

4  52136 

5  542136 

6  521436 
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Figure  7:  Illustration  of  routing  for  telephone  network  example 


fn 

fl2 


/lR(l) 


fmR(m)  ’ 

51 

52 


c2  := 

Q(m7^(m))xl 

jmxl 

,  b2:= 

e 

a 

QTl  X  1 

a 

Sm 

^2 


ail  ai2 

•  ■  •  a^i) 

QTiXm 

jnxn 

An  *  = 

-^lxfi(l) 

QlX^i) 

q1  x  R(2)  ... 

^lxR(2) 

^*•2  ■  — 

0lx-R(1) 

Qlx«(2)  ... 

-^lxtt(m) 

jmxm 

Qmxn 

To:-- 


_ jnxn  Q 


0" 


0 


qc; 


12  Bond  investment  planning 

Due  to  K.  Frauendorfer,  C.  Marohn,  M.  Schurle  [6] 
(Multistage,  linear  stochastic  problem) 


12.1  Description 

Frauendorfer,  Marohn,  and  Schurle  [6]  describe  a  suite  of  test  problems  for 
multistage  stochastic  programming,  based  on  bond  investments.  The  test 
problems  are  denoted  SGPFmYn,  where  m  E  {3, 5},  and  n  E  {3, 4, 5, 6, 7}. 

Many  business  ventures  are  financed  by  lending  bonds,  and  many  of  these 
ventures  also  purchase  bonds.  There  is  risk  in  such  dealings,  as  returns  on 
bonds  fluctuate,  and  earnings  from  the  business  ventures  are  uncertain.  This 
risk  cannot  be  modeled  by  assuming  a  mean  rate  of  return.  Therefore,  the 
scenario  is  a  good  one  for  the  application  of  stochastic  programming. 

Bonds  mature  in  certain,  standard  time  periods.  Suppose  we  will  consider 
transactions  in  bonds  with  standard  maturities  in  the  set  Vs.  Suppose  that 
the  longest  maturity  in  Vs  is  D  months.  Then,  since  the  time  frame  is  rolling 
in  such  problems,  we  must  include  in  the  model  bonds  which  mature  in  d 
months,  where  d  EV  =  { 1,2,...  ,  D}. 

Let  vf’+  be  the  amount  of  new  borrowing  done  at  time  t  with  maturity 
d  E  Vs,  and  let  vf~  be  the  amount  of  new  lending  done  in  the  same  cir¬ 
cumstances.  Then,  if  vf  is  the  balance  of  bond  transactions  at  time  t  with 
maturity  d,  we  have 

vf’~  if  d  E  Vs , 

if  dEV\Vs. 

The  total  balance  of  bond  transactions  at  time  t  is 

xt  =  J2vt- 

dev 

If  this  quantity  is  positive,  the  balance  will  be  used  to  fund  the  business 
venture  during  the  time  period  t.  Rather  than  writing  xt  as  a  function  of 
historical  balances  and  rates  of  return,  Frauendorfer,  Marohn,  and  Schurle 
[6]  simply  express  it  as  the  stochastic  quantity 

Xt  =  Xt- 1  +  £t, 


vt  = 


vts  +  vf'* 

ft1 
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where  £t  is  a  random  variable. 

To  limit  the  sale  of  bonds,  the  authors  [6]  include  the  constraint 


devs 


X>-i 


d=  1 


<&. 


Given  random  rates  of  return  rjd,+,  and  r]*,  corresponding  to  the 
quantities  vf~,  vf’+,  and  xt,  respectively,  the  objective  is  to  maximize  the 
expected  return: 


E 


-rii’+v{>+}+r£xt 


The  returns  at  time  t  =  0  are  actually  deterministic.  So  the  decision  variables 
for  time  t  —  0  are  the  so-called  first  stage  decision  variables  in  the  stochastic 
problem. 


12.2  Problem  statement 

We  change  the  problem  to  a  minimization,  and  separate  the  first  stage  vari¬ 
ables  and  constraints  from  the  recourse  variables  and  constraints.  We  are 
given  all  values  for  the  time  t  =  —  1  decision  variables,  and  the  time  t  =  0 
values  of  all  77  and  £.  The  program  then  is  to 

minimize  ^  [— 7?o’  +  Vo’+vo+]  ~  + 

devs 

E  [  E  +  vt’+vt’+]  -  r%Xt 

t=i  \<zex>s 

subject  to 

vd  -  ulj1  -  v%'+  +  Vq’~  =  0,  We  Vs,  (41) 

vd  -  v 1+1  =0,  We  V\DS,  (42) 

dev 

Xq  37—1  Co, 


97 


(43) 


M 

e  ”o'+ -  i>->  -  &>’ 

devs  d=i 

v%'+X'~  >  0,  WeV, 

\/devs,t  =  l,...  ,T, 

VdeV\Ds,t  =  1,...  ,T, 

Vt  =  l,...,T, 

Vi  =  1, . . .  ,T, 

Vi  =  1, . . .  ,T,  (44) 

\/d  £V,t  —  1, . . .  ,T. 


12.3  Numerical  examples 

A  total  of  ten  numerical  examples  in  SMPS  format  [1]  are  available  from 
Birge’s  POSTS  web  site  [9].  Since  the  only  coefficients  to  be  specified  in 
this  model  are  stochastic,  specifying  any  one  problem  here  would  require 
duplicating  the  stochastic  file  from  the  set  of  SMPS  files.  Therefore,  we  refer 
the  reader  to  the  publicly  available  SMPS  files  [9]. 

12.4  Notational  reconciliation 

We  may  rearrange  the  equations  represented  by  (41)  and  (42)  so  that  they 
are  in  ascending  order,  by  d.  Then  we  have  D  constraints,  each  with  right- 
hand  sides  v^1,  and  left  hand  sides  depending  on  whether  d  G  Vs  or  not. 
We  replace  all  vf  with  the  term  ( vp f  —  vmf),  with  the  added  constraints  that 
vpf,  vmf  >  0. 

Let  {dl,  d2, . . .  ,  dN}  -  {d  G  Vs}.  We  define  the  matrix  A5  G  RDxN  by 
As  :=  [  edl  ed 2  •  •  •  edN  ]  , 
where  e1  G  M.D  is  the  ith  unit  vector. 


vf  -  vd+l  -  vf’+  +  vf'  =  0, 


n  ~  vt-i 


d+ 1 


U. 


xt  -  Yv*  =  °> 

dev 

%t  %t— i 
M 

dev s  d—  1 

vdt'+yt'-  >  o, 
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Let  st  be  the  slack  variable  associated  with  constraint  (43)  or  (44).  Assign 
the  new  notation 


by  := 


(vp2_i  —  vm2_t) 

(vpz_x  —  vm lx) 

(vpEx  -  wm®) 

0 
0 

6  +  x-i 

to  +  ZhivpU-vrnU)  J 


Ci  := 


-Po 

0 


4N’+ 


~VoN’ 

q2DxO 


and 

aa  :=  [  «o 
Also,  let 


dl,-f 


dNt+  nd\- 


dNr 


vp\  •• 

•  VPo 

vm\  •  •  ■ 

-  0Dxl 

0Dxl 

— As 

As 

jDxD 

_ jDxD 

1 

0 

0!XN 

olxN 

_jlx£) 

1  lxD 

A\ 

1 

0 

q1  x  N 

glxiV 

o lxD 

olxD 

0 

1 

^IxN 

O 

x 

0lxD 

0lxD 

vmt 


?r 


Analogous  assignments  are  made  for  t  =  2,3, ...  ,T,  except  that  ct  is 
made  stochastic  for  these  times,  because  rj  is  stochastic.  Also, 


bt  := 


0Dxl 

0 

0 

6-1 

6-i 


and 


Tt  := 


|£»xl 

o 

t) 

x 

q  DxN 

qDxN 

_ jDxD 

jDxD 

0 

0 

qIxN 

o 

>— 1 

X 

olxD 

o lxD 

-1 

0 

0lxN 

o 

X 

QiXO 

X 

y—{ 

O 

0 

0 

q1  xN 

qI  xN 

-Wi 

Wi 
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where 


Wi  :=  [  llxM  olx^_M^  ]  . 
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